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Executive  Summary 


This  is  the  Final  Report  on  the  AFOSR  project  (FA9550-11-1-0111)  entitled: 

Physics  based  modeling  of  compressible  turbulence.  The  period  of  performance  was,  June 
15,  2011-June  14,2016. 

The  following  Postdoctoral  Fellows  and  students  of  the  Center  for  Turbulence  Research 
(CTR)  contributed  to  this  project: 

Dr.  Johan  Larsson  (currently,  Assistant  Professor,  University  of  Maryland) 

Dr.  Joseph  Nichols  (currently,  Assistant  Professor,  University  of  Minnesota) 

Dr.  Ivan  Bermejo-Moreno  (currently,  Assistant  Professor,  University  of  Southern  California) 
Dr.  Taraneh  Sayadi  (currently,  Assistant  Professor,  University  of  Illinois, 
Urbana/Champaign) 

Dr.  Ik  Jang,  Stanford  University  (Ph.D.,  August,  2016). 

Curtis  Hamman,  Stanford  University 

Important  challenges  in  high  fidelity,  high-speed  flow  simulations  were  considered  in  this 
project,  including:  shock/turbulence  interaction,  aerodynamic  noise,  supersonic  propulsion 
and  boundary  layer  transition.  The  latter  is  being  continued  and  expanded  with  support 
from  a  new  AFOSR  contract.  This  report  focuses  on  the  first  three  challenges. 

The  following  are  the  project’s  key  accomplishments: 

•  Validation  of  wall  models  in  numerical  simulation  of  shock/boundary  layer 
interaction 

•  Demonstration  of  sidewall  confinement  effects  in  shock/boundary  layer  interaction 

•  LES  validation  of  jet  noise  reduction  with  chevrons 

•  Mechanism  of  crackle  in  supersonic  jet  noise 

•  Direct  numerical  simulation  of  H-type  and  K-type  transition  to  turbulence 

•  The  first  LES  computation  using  one  million  computer  cores 

•  Stability  analysis  of  scramjet  unstart 

•  Wall  modeled  LES  simulation  and  validation  of  HIFIRE  scramjet  engine 
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Part  1 

Physics-based  modeling  of  compressible  turbulence 

1.1  Validation  of  wall  models  in  large-eddy  simulation  of  shock¬ 
boundary  layer  interactions  with  sidewall  confinement 

Despite  the  rich  history  of  wall-modeling  developments  in  LES  found  in  the  literature,  its 
application  to  cases  with  shock  waves  is  scarce.  The  study  presented  in  this  section  stems 
from  the  recent  experiments  ofHelmer  etal.  [2012]  and  Campo  etal.  (2012]  and  aims  to 
validate  and  use  wall-modeled  large-eddy  simulations  (WMLES)  to  explore  the  physics  of  a 
series  of  shock/turbulent-boundary-layer  interactions  (STBLIs)  of  increasing  strength 
inside  a  duct  with  a  low  aspect  ratio  cross-section.  The  confinement  effects  imposed  by  the 
side  walls  in  a  low-aspect  ratio  duct  configuration  are  also  studied,  by  comparing 
simulations  that  include  those  side  walls  with  spanwise  periodic  simulations. 

Despite  all  the  previous  work  devoted  to  the  study  of  the  low-frequency  unsteadiness  in 
STBLI,  little  is  known  about  the  influence  (if  any]  of  the  confinement  imposed  by  the 
sidewalls.  Dussauge  etal.  (2006)  considered  several  experiments  of  reflected  STBLIs  and 
hypothesize  that  the  three-dimensional  structure  of  the  separation  bubble  might  be  the 
cause  of  the  low-frequency  unsteadiness.  Another  objective  of  this  work  is  thus  to  study  the 
effect  of  sidewalls  on  STBLI  low-frequency  motions. 

The  results  presented  in  this  section  are  further  detailed  in  Bermejo-Moreno  et  al.  (2014). 

1.1.1  Computational  setup  and  methodology 

The  computational  domain  is  shown  in  Fig.  1.  It  includes  the  test  section  of  the  wind  tunnel 
used  in  the  experiments  byHelmer  &  Eaton  (2011),  Helmer  etal.  (2012)  and  Campo 
et  al.  (2012),  consisting  of  a  constant-area  duct  section  (45x47.5  mm2)  followed  by  a  short 
contraction  produced  by  a  compression  wedge  that  spans  the  top  wall  at  a  20°  angle,  and 
another  constant-area  section  downstream  of  the  apex  of  the  compression  wedge.  The 
reference  system  is  chosen  such  that  the  origin  of  the  streamwise  coordinate,  x,  coincides 
with  the  location  of  the  foot  of  the  compression  wedge;  y  is  the  vertical  coordinate,  with 
origin  on  the  bottom  wall;  and  z  is  the  spanwise  coordinate,  with  origin  on  the  left  wall  as 
one  looks  downstream  from  inside  the  duct.  The  compression  wedge  (nominally  two- 
dimensional)  interacts  with  the  incoming  boundary  layer  developed  on  the  top  wall, 
generating  an  oblique  incident  shock  wave  that  reflects  off  the  bottom  wall  boundary  layer. 
In  the  present  study,  we  consider  independently  three  different  heights  of  the  compression 
wedge:  1.1,  3  and  5  mm.  The  strength  of  the  incident  shock  and  of  the  subsequent 
interactions  with  the  boundary  layers  increases  with  the  height  of  the  compression  wedge. 
At  the  apex  of  the  compression  wedge,  an  expansion  fan  is  generated,  turning  the  flow  back 
to  follow  the  horizontal  top  wall.  This  expansion  fan  propagates  downstream  and 
eventually  interacts  with  the  flow  along  the  bottom  wall. 
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Fig.  1.  Computational  domain  and  schematic  representation  of  incident  (ISW)  and  reflected 
(RSW)  shock  waves  and  turbulent  boundary  layers  (TBL). 


The  present  simulations  solve  the  spatially-filtered  compressible  Navier-Stokes  equations 
for  the  conserved  variables  of  mass,  momentum  and  total  energy  of  a  calorically  perfect 
gas,  using  a  finite  volume  formulation,  control-volume  based  discretization  on 
unstructured  hexahedral  meshes.  A  solution-adaptive  methodology  is  implemented,  that 
combines  a  non-dissipative  centered  numerical  scheme  and  an  essentially  non-oscillatory 
(ENO)  second-order  shock-capturing  scheme.  The  latter,  which  uses  an  HLLC  Riemann 
solver  for  the  computation  of  Euler  fluxes,  is  applied  in  regions  near  shock  waves, 
identified  by  a  shock  sensor  activated  according  to  the  criterion:  -dui*0Xk  >  max(^^, 
0.05oA),  where  dui^dxk,  cojOOj  and  c  are  the  local  dilatation,  enstrophy  and  sound  speed, 
respectively,  and  A  is  the  mesh  cell  size.  A  mesh-based  blend  of  the  centered  and  upwind 
numerical  schemes  is  used  for  numerical  robustness  (see  Khalighi  etal.,  2011).  Subgrid 
scale  stresses  are  explicitly  modeled  following  Vreman  (2004],  and  a  fixed  turbulent 
Prandtl  number  of  0.9  is  used  to  model  the  subgrid  scale  heat  flux.  For  each  wedge  height, 
three  independent  simulations  were  run  on  meshes  of  increasing  resolution  to  perform  a 
grid-convergence  study.  The  wall  model  proposed  by  Kawai  &  Larsson  (2012)  is  used.  It 
solves  the  equilibrium-boundary-layer  equations  in  a  refined,  near-wall  inner  grid, 
embedded  in  the  coarser,  background  LES  grid  The  wall  model  is  applied  at  all  four  walls, 
considered  adiabatic. 

To  account  for  the  turbulent  nature  of  the  boundary  layers  at  the  inlet  of  the  simulations, 
we  use  a  synthetic  turbulent  inflow  generator  based  on  a  digital  filtering  technique 
originally  proposed  by  Klein  etal.  (2003),  with  the  improvements  ofXie  &  Castro  (2008) 
and  Touber  &  Sandham  (2009).  The  availability  of  experimental  PIV  data  on  three  vertical 
planes  near  the  sidewall,  in  addition  to  the  near-center  plane,  makes  the  characterization  of 
the  inflow  more  complete  than  what  is  commonly  available  in  previous  STBLI  experiments. 
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1.1.2  Results 


1.1.2.1  Validation  with  experimental  data 

To  compare  with  experimental  PIV  measurements,  flow  variables  in  the  WMLES  were 
averaged  in  time  over  a  period  of  10005(/iJo  after  the  initial  transient.  The  reflected  STBLI 
region  near  the  bottom  wall  on  the  spanwise-normal  center  plane  of  the  duct  is  enlarged  in 
Fig.  2,  for  both  wedge  heights.  Color  contours  of  mean  streamwise  and  vertical  velocity 
extracted  from  the  PIV  are  shown  in  the  background,  superimposing  the  corresponding 
contour  lines  from  the  WMLES.  A  reasonable  agreement  is  observed  between  experiments 
and  simulations  for  these  mean  quantities  in  the  plane  under  consideration. 

A  representative  quantitative  comparison  of  PIV  and  WMLES  mean  velocity  profiles  for  the 
3  mm  compression  wedge  is  shown  in  Fig.  3  for  the  three  mesh  resolutions  considered,  so 
that  the  grid-convergence  of  the  results  can  be  assessed.  The  shaded  grey  area  in  each  plot 
represents  the  wall-normal  extent  of  the  inner  grid  used  in  the  simulation  by  the  wall- 
model.  We  note  that  the  simulation  profiles  shown  in  these  figures  correspond  to  the  LES 
background  grid  solution.  Inside  the  grey  shaded  regions,  the  wall-model  solution  [not 
shown)  is  being  used  instead  by  the  simulation  to  provide  the  wall  boundary  condition  to 
the  LES.  Similar  qualitative  and  quantitative  agreement  between  experiments  and 
simulations  is  found  for  turbulence  quantities  (not  presented  in  this  report). 


Fig.  2.  Contours  of  mean  streamwise  (top)  and  vertical  (bottom)  velocity  on  the  first 
reflected  STBLI  off  the  bottom  wall,  taken  on  a  vertical  xy-plane  located  at  7/ 60  =  3.89  (near 
the  center  of  the  duct)  for  the  1.1  mm-high  wedge  case  (left)  and  3  mm-high  wedge  case 
(right).  PIV  data  shown  in  the  background  (by  a  colour  map)  with  24  contour  lines  from  the 
WMLES  data  superimposed  (solid  and  dashed  lines  correspond  to  positive  and  negative 

isocontour  values,  respectively). 
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(xicr1)  (xlO-1)  (xio-1)  (xlO-1)  (xlO-1)  (xlO-1)  (xlO-1)  (xlO-1)  (xlO-1)  (xlO-1) 

Fig.  3.  Mean  velocity  profiles  in  the  reflected  STBLI  for  the  3  mm-high  wedge  case:  (a]  and 
(b),  streamwise  velocity  (IRJ0);  (c)  and  (d),  vertical  velocity  (V4J0).  [a)  and  (c):  near-center 
plane  (2/ S0  =  3.89];  (b)  and  (d):  near-sidewall  plane  (z4>0  ~  0.7).  From  left  to  right  in  each 
plot:  x/50  =  2.1,  4.6,  5.4,  6.4,  7.9.  Symbols:  PIV  data;  solid,  WMLES  fine  resolution;  dashed, 
WMLES  medium  resolution;  dotted,  WMLES  coarse  resolution.  Shaded  grey  areas  near 
y/60  =  0  represent  the  extent  of  the  inner,  wall-model  layer  in  the  simulations. 

1.1. 2. 2  Confinement  effects  on  the  strong  interaction 

Once  confidence  in  the  simulation  methodology  was  built  for  the  1.1  mm-  and  3  mm-high 
wedge  cases,  additional  simulations  were  performed  for  a  stronger  interaction  case 
corresponding  to  a  5  mm-high  wedge  configuration,  keeping  the  same  20°  wedge  angle. 
The  primary  goal  of  these  simulations  was  to  produce  a  larger  region  of  mean  flow  reversal 
by  an  increased  shock  strength,  to  study  the  effects  of  the  confinement  imposed  by  the 
sidewalls  on  the  extent  and  dynamics  of  the  separation  bubble. 

Two  sets  of  simulations  with  the  5  mm-high  wedge  were  performed:  one  set  includes  the 
sidewalls  whereas  the  other  set  uses  spanwise  periodic  boundary  conditions,  resulting  in  a 
nominally  two-dimensional  interaction.  Three  meshes  of  increasing  resolution  were 
considered  for  each  set.  The  streamwise  length  of  the  computational  domain  in  both  cases 
is  identical  to  the  3  mm-high  wedge  configuration.  The  spanwise-periodic  simulations  have 
a  reduced  domain  width  of  4.460  with  a  uniform  grid  spacing  in  z  (i.e.,  no  stretching  away 
from  the  center  plane). 

Fig.  4  shows  contours  of  streamwise  and  vertical  mean  velocities  for  simulations  with  and 
without  sidewalls.  Mean  sonic  and  separation  lines  are  superimposed  as  dashed  and  solid 
lines,  respectively,  in  each  plot.  Regions  of  mean  flow  separation  appear  in  both 
simulations  in  the  compression  STBLI  [top  wall)  and,  as  targeted,  in  the  first  reflected 
STBLI  (bottom  wall).  The  most  striking  difference  between  the  two  simulations  is  the 
character  of  the  first  reflected  STBLI  at  the  bottom  wall:  the  spanwise-periodic  simulation 
shows  a  regular  shock  intersection  whereas  the  simulation  with  sidewalls  presents  a 
singular  shock  intersection  centered  around  (x,y)  «  (3.43,  2.75)60.  In  the  singular 
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intersection,  a  Mach  stem  is  formed  in  between  the  incident  and  reflected  shock  waves, 
consisting  of  a  quasi-normal  shock  with  two  triple  points  at  the  intersection  with  the 
incident  and  separation  shocks  and  the  reflected  and  transmitted  shocks, 
respectively  (Babinsky  &  Harvey,  2011].  A  region  of  subsonic  flow  appears  immediately 
downstream  of  the  Mach  stem,  which  is  accelerated  by  the  surrounded  regions  of 
supersonic  flow  on  each  side  until  the  flow  in  the  wake  of  the  Mach  stem  becomes 
supersonic  again,  passing  through  a  sonic  throat.  In  the  spanwise-periodic  simulation  the 
Mach  stem  is  not  present  and  a  regular  shock  intersection  occurs,  in  contrast  with  the 
simulation  with  sidewalls.  This  indicates  that,  for  this  configuration  and  flow  parameters, 
confinement  effects  imposed  by  the  sidewalls  are  responsible  for  strengthening  the 
incident  and  separation  shocks  beyond  which  their  two  polars  cannot  intersect,  requiring  a 
quasi-normal  shock  and  two  triple  points  that  provide  the  singular  Mach  intersection. 

Both  spanwise-periodic  and  sidewall  simulations  present  regions  of  mean  flow  reversal  in 
the  compression  wedge  interaction  at  the  top  wall  as  well  as  in  the  first  reflected  STBLI  at 
the  bottom  wall.  Fig.  5  (a)  shows  a  three-dimensional  visualization  of  the  separation  bubble 
on  the  bottom  wall.  The  Mach  stem  and  the  structure  of  the  transmitted  and  reflected 
shocks  of  that  STBLI  are  also  shown,  along  with  the  bottom  wall  (colored  by  mean 
pressure)  and  the  vertical  plane  at  the  center  of  the  duct  (z^o  =  4.4),  colored  by  mean 
vertical  velocity,  for  reference.  Fig.  5(b)  shows  the  projection  of  several  slices  of  the 
separation  bubble  on  vertical  planes  normal  to  the  spanwise  coordinate,  at  increasing 
distance  from  the  center  of  the  duct.  A  dashed  line  corresponding  to  the  mean  separation 
bubble  for  the  spanwise-periodic  simulation  is  also  shown  for  comparison. 


0.0  0.2  0.5  0.8  1.0  -0.2  -0.1  0.0  0.1  0.2 

urn.  vm. 

Fig.  4.  Contours  of  mean  streamwise  (a,c)  and  vertical  (right)  velocities  for  the  5  mm-high 
wedge  case  at  the  center  plane  (z4)o  =  4.4),  with  sonic  line  (dashed)  and  separation  lines 
(solid)  superimposed  in  each  plot.  Comparison  between  simulations  with  sidewalls  (top) 
and  spanwise  periodicity  (bottom).  I:  incident  shock  generated  by  compression  corner;  MS: 
Mach  stem;  S,  T  and  R:  separation,  transmitted  and  reflected  shocks  from  first  STBLI  at  the 
bottom  wall;  R':  reflected  shock  from  the  second  STBLI  at  the  top  wall;  Q:  incident  shock  on 
the  third  STBLI  at  the  top  wall;  R reflected  shock  from  second  STBLI  at  bottom  wall. 
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Fig.  5.  [a)  3D  visualization  near  the  first  STBLI  at  the  bottom  wall  for  the  5  mm-high  wedge 
simulation,  showing  the  separation  bubble  (isosurface  of  zero  mean  streamwise  velocity 
colored  by  vertical  Reynolds  stress,  vVv'/U0  from  0  to  0.1,  black  to  white),  the  Mach  stem 
with  the  reflected  shocks  educed  by  an  isosurface  of  constant  mean  density  p/p0  =  2.4, 
colored  by  mean  temperature  from  TT0  =  1.37  to  1.53  (white  to  black),  the  center  vertical 
plane  (z4>0  =  4.4)  with  isocontours  of  vertical  mean  velocity,  V4J0,  from  -0.2  to  0.2  (blue  to 
green  to  red,  semitransparent),  and  the  bottom  wall  (with  the  full  spanwise  extent  between 
sidewalls,  colored  by  mean  pressure,  pt>0  from  1.0  to  2.1,  blue  to  gray  to  red),  (b)  Isolines  of 
zero  mean  streamwise  velocity  (defining  the  mean  separation  bubble)  for  the  simulation 
with  sidewalls  (solid)  on  z-normal  planes  from  the  center  of  the  duct  towards  the  sidewall 
(in  the  direction  of  the  arrow),  and  for  the  spanwise-periodic  simulation  (dashed  line). 

The  separation  bubble  for  the  spanwise-periodic  simulation  is  located  approximately  50 
farther  downstream  and  is  smaller  (about  20%  shorter  streamwise  and  of  nearly  half  the 
height)  than  in  the  simulation  with  sidewalls,  whereas  the  skewed  shape  is  similar  in  both 
simulations.  The  smaller  size  found  for  the  spanwise-periodic  simulation  is  consistent  with 
the  observations  byPriebe  etal.  (2009),  Hadjadj  etal.  (2010),  Pirozzoli  etal.  (2010), 
Pirozzoli  &  Bernardini  (2011)  and  Morgan  etal.  (2013),  who  found  an  underprediction  of 
the  streamwise  separation  length  when  comparing  spanwise-periodic  simulation  results 
with  experiments. 

Fig.  6  shows  streamwise  profiles  of  the  mean  streamwise  skin  friction  coefficient,  Cf,  and 
the  wall  pressure,  pw,  obtained  on  the  bottom  wall  at  several  planes  normal  to  the  spanwise 
direction.  Profiles  for  the  spanwise-periodic  simulation  are  also  plotted  in  dashed-dotted 
lines.  The  mean  separation  bubble  shows  a  strong  three-dimensionality  imposed  by  the 
lateral  confinement.  For  the  spanwise-periodic  simulation,  the  separation  bubble  is 
significantly  smaller,  in  agreement  with  earlier  findings  in  the  literature.  The  profiles  of 
streamwise  skin  friction  coefficient  and  wall  pressure  for  the  interaction  strengths 
simulated  confirm  a  link  between  the  three-inflection-point  wall-pressure  profile  and 
presence  of  mean  flow  reversal. 
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Fig.  6.  Profiles  of  time-averaged  streamwise  skin  friction  coefficient  Cf  [a),  and  wall 
pressure  pw  (b),  along  the  bottom  wall,  for  the  5  mm-high  wedge  simulation  with  sidewalls: 
solid  line,  center  plane  (#80  =  4.4);  dashed  lines  of  decreasing  dash  length  correspond  to 
decreasing  distance  to  the  sidewall  (z4>0  =  3.85,  3.32,  2.79,  2.25, 1.72, 1.19,  0.66,  0.13).  The 
spanwise-averaged  profiles  retrieved  from  the  spanwise-periodic  simulation  are  plotted  in 

dash-dotted  lines  for  comparison. 


1. 1.2.3  Low-frequency  unsteadiness 

Fig.  7  shows  contours  of  the  premultiplied  power  spectral  density  (PSD  )  of  the  bottom-wall 
pressure  signals  as  a  function  of  the  streamwise  location  and  the  normalized  frequency 
f6o4J0.  The  PSD'  at  each  streamwise  location  is  normalized  with  its  integrated  value  across 
all  frequencies.  Fig.  15(a)  and  (b)  correspond  to  the  simulation  with  sidewalls,  at  spanwise 
locations  of  z4)0  =  4.4  (center  plane)  and  t/50  ~  0.5  (near-sidewall  plane),  respectively. 
Fig.  7(c)  shows  the  centerplane  results  for  the  spanwise-periodic  simulation. 

Near  the  first  STBLI  on  the  bottom  wall  (x4>0  *  2),  regions  of  low-frequency  content  are 
visible  in  the  wall-pressure  spectra  in  Fig.  7(a,b,c),  extending  to  frequencies  below  St  «  0.1 
and  peaking  around  St  ~  0(0.01).  At  the  center  plane  of  both  simulations  (Fig.  7a  and  c), 
the  streamwise  location  where  these  low-frequency  signals  are  found  coincides  with  the 
foot  of  the  separation  shock,  located  upstream  of  the  mean  separation  point.  The  spanwise- 
periodic  simulation  shows  a  downstream  shift  of  the  location  with  low-frequency  signals  of 
approximately  60,  consistent  with  the  delayed  position  of  the  first  STBLI  on  that  bottom 
wall  compared  with  the  simulation  with  sidewalls. 

Moving  downstream  along  the  first  STBLI,  a  shift  is  observed  in  the  spectrum  from  the 
characteristic  frequencies  of  the  incoming  turbulent  boundary  layer  toward  lower 
frequencies.  The  spanwise-periodic  simulation  shows  a  similar  behavior  inside  the  first 
STBLI  region  than  already  described  for  the  center  plane  of  the  simulation  with  sidewalls, 
the  only  appreciable  difference  being  the  reduced  streamwise  extent  of  the  region  of  high 
intensity  spectral  density,  as  it  begins  farther  downstream  (xtf>0  ~  3)  than  for  the  simulation 
with  sidewalls  (x4>0  *  1.5),  ending  approximately  at  the  same  location  in  both  cases  (xtf>0~8). 
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Low-frequency  unsteadiness  reappears  at  the  second  STBLI  on  the  bottom  wall  (see  Fig.  7), 
particularly  at  the  center  plane  (xtf>0 »  13  for  the  simulation  with  sidewalls  and  #S0  ~  15  for 
the  spanwise  simulation].  The  spectral  density  content  at  this  second  STBLI  is  dominated 
by  such  low  frequencies  (St  G  0.01  -  0.1),  suggesting  a  much  stronger  unsteadiness  than  for 
the  first  STBLI.  Near  the  sidewall,  a  trace  of  the  low-frequency  unsteadiness  derived  from 
the  second  STBLI  is  captured  in  the  PSD'  at  #S0  ~  13,  although  its  intensity  is  much  weaker 
than  at  the  center  plane. 

1.2  Wall  modeled  LES  and  validation  of  HIFiRE  scramjet 

As  part  of  our  validation  plan  for  the  physics-based  modeling  of  compressible  flows  we 
conducted  large-eddy  simulation  of  the  HIFiRE-2  scramjet  at  two  operating  conditions 
using  an  equilibrium  wall  model  and  flamelet-based  combustion  models. 

The  Hypersonic  International  Flight  Research  Experimentation  Program  (HIFiRE)  is  a 
collaborative  effort  among  the  United  States  Air  Force  Research  Laboratory  (AFRL),  NASA, 
and  the  Australian  DSTO.  Its  second  flight  test  (HIFiRE-2)  was  successfully  flown  in  May 
2012.  Key  elements  that  set  the  HIFiRE-2  scramjet  apart  from  previous  scramjet 
experimentation  programs,  posing  additional  simulation  challenges,  are:  1)  mode 
transition,  that  is,  the  operation  at  variable  Mach  number  from  dual  (ramjet)  mode  to 
scramjet  mode  operation  over  a  flight  Mach  range  from  6  to  8  at  nearly  constant  dynamic 
pressure;  2)  the  use  of  a  hydrocarbon  fuel,  instead  of  hydrogen;  3)  a  multi-staged  fuel 
injection  system;  and  4)  the  presence  of  a  cavity-based  flame-holder  located  between  the 
two  injection  stages. 

HIFiRE-2  was  supported  by  a  campaign  of  ground  tests  (Hass  et  al.,  2009;  Cabell 
et  al.,  2011)  performed  in  the  HIFiRE  Direct  Connect  Rig  (HDCR)  at  the  NASA  Langley  Arc- 
Heated  Scramjet  Test  Facility  (AHSFT).  These  ground  tests,  for  which  public  data  are 
available  for  both  the  dual  and  scramjet  modes  of  operation,  are  the  target  of  our 
simulations.  In  addition  to  the  subgrid  scale  and  wall  models  already  considered  in  the 
context  of  LES  of  STBLI,  these  simulations  include  combustion  models.  The  complexity  of 
the  chemical  mechanism  of  a  hydrocarbon  fuel  such  as  the  one  used  in  HIFiRE-2,  involving 
a  large  number  of  species  and  chemical  reactions,  also  adds  to  the  need  for  a  tractable 
combustion  modeling  approach. 


Further  details  of  the  study  reported  in  §2.2  can  be  found  in  Bermejo-Moreno  et  al.  (2013). 
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Fig.  7.  Contour  maps  of  premultiplied  power  spectral  density  (PSD']  (arbitrary  scale)  of 
bottom-wall  pressure  for  an  array  of  probes  located  along  the  streamwise  coordinate:  (a) 
center  plane  ( t/50  =  4.4)  of  simulation  with  sidewalls,  (b)  near-sidewall  plane  (z4>0  *  0.5)  of 
simulation  with  sidewalls,  (c)  center  plane  of  spanwise-periodic  simulation.  PSD's  are 
normalized  for  each  streamwise  coordinate  location. 
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1.2.1  Experimental  geometry  and  flow  conditions 

Fig.  8  shows  a  side  view  of  the  geometry  of  the  HIFiRE-2  scramjet  engine,  including  the 
isolator  and  combustor  stages.  The  overall  length  is  711.3  mm,  with  a  constant  width  of 
101.6  mm.  The  isolator  comprises  nearly  one  third  of  the  total  length  of  the  engine  and  has 
a  constant  height  of  25.4  mm,  whereas  the  top  and  bottom  walls  [body  and  cowl  sides, 
respectively)  of  the  combustor  diverge  at  a  constant,  total  angle  of  2.6°.  Two  opposed  cavity 
flameholders  are  located  in  the  combustor,  separating  the  two  fuel  injection  stages,  each 
with  eight  injection  ports  [four  on  the  body  side  and  four  on  the  cowl  side,  equispaced 
spanwise).  The  primary  injectors,  located  upstream  of  the  cavity,  are  angled  at  a  15° 
inclination  from  the  wall,  with  a  circular  diameter  of  3.175  mm.  Secondary  injectors, 
located  downstream  of  the  cavity,  are  perpendicular  to  the  wall  with  a  2.38  mm  circular 
diameter. 


Isolator  (203.2  mm) 

f - — - f 


Combustor  (508.1  mm) 


1=7- 

I - ►X 


Primary  injectors 


1 


T 


Cavity  Secondary  injectors  1.3°  expansion 

Fig.  8.  Schematic  representation  of  the  HIFiRE  scramjet,  including  the  isolator  and 

combustor.  Flow  from  left  to  right. 


Two  experimental  conditions  are  considered  in  our  simulations,  as  summarized  in  Table  1. 
These  conditions  correspond  to  the  experimental  runs  labeled  125.1  and  136.3  (see  Hass 
et  al.,  2009;  Cabell  etal.,  2011)  for  the  simulated  flight  Mach  numbers  of  6.5  and  8.0, 
respectively,  representative  of  the  two  operational  modes  of  the  HIFiRE-2  system  (i.e., 
dual-mode  and  scram-mode). 


Test  ID  Mf  Mn  mair  rhfi  m/2 


125.1  6.5  2.51  1.137  0.0291  0.0437 

136.1  8.0  3.46  0.912  0.0234  0.0354 

Table  1.  Flow  configurations  for  the  two  ground  tests  targeted  in  the  simulations.  Mf, 
simulated  flight  Mach  number;  Mn,  HDCR  nozzle  Mach  number;  mair,  air  mass  flow  rate;  mn 
and  ihf2,  fuel  mass  flow  rate  in  the  primary  and  secondary  injection  stages,  respectively. 

Mass  flow  rates  in  kg/s. 


1.2.2  Modeling  approach 

The  spatially  filtered  compressible  Navier-Stokes  equations  for  the  conserved  variables  of 
mass,  momentum,  and  total  energy  (which  includes  sensible,  kinetic,  and  chemical  energy 
for  the  reacting  flow)  are  solved  in  the  present  simulations  using  a  finite  volume 
formulation  on  unstructured  hexahedral  meshes.  The  solver  implements  a  solution- 
adaptive  methodology  that  combines  a  non-dissipative  centered  numerical  scheme  and  an 
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essentially  non-oscillatory  (ENO)  second-order  shock-capturing  scheme  [with  an  HLLC 
Riemann  solver  for  the  computation  of  the  Euler  fluxes],  which  is  applied  only  in  regions 
near  shock  waves  identified  by  a  shock  sensor  activated  according  to  the  criterion:  -dui^dxk 
>  rnaxf^^,  0.05oA),  where  dui*0Xk,  (UjOOj  and  c  are  the  local  dilatation,  enstrophy,  and 
sound  speed,  respectively,  and  A  is  the  control  volume  size.  Additionally,  the  shock¬ 
capturing  scheme  is  applied  in  cells  where  differences  with  adjacent  cells  in  mixture 
fraction  or  temperature  are  greater  than  0.4  and  2500  K,  respectively.  Away  from 
discontinuities,  the  non-dissipative  [second-order]  scheme  is  applied.  A  mesh-based  blend 
of  centered  and  upwind  numerical  schemes  is  used  for  robustness  [Khalighi  etal.,  2011]. 
Subgrid  scale  stresses  are  modeled  following  Vreman  [2004].  Gradient-diffusion  models 
are  used  for  the  subgrid  scale  heat  flux  and  species  transport  with  fixed  turbulent  Prandtl 
and  Schmidt  numbers  of  0.9.  Subgrid-scale  model  terms  are  set  to  zero  where  the  shock¬ 
capturing  scheme  is  active  [i.e.,  in  regions  marked  by  the  shock  sensor],  to  avoid  adding 
extra  dissipation  to  the  already  dissipative  ENO  scheme  [Bermejo-Moreno  etal.,  2010].  A 
three-stage,  third-order  explicit  Runge-Kutta  algorithm  is  used  to  advance  the  discretized 
equations  in  time. 

Meshes  of  increasing  resolution  up  to  700  million  control  volumes  were  used,  with  grid 
spacing  refined  near  the  walls,  fuel  injectors  [meshed  using  O-grids],  and  in  the  region  of 
development  of  the  shear  layer  in  the  cavity. 

The  wall  model  proposed  byKawai  &  Larsson  [2012]  is  also  used  in  the  present 
simulations.  As  stated  in  the  previous  section,  this  wall  model  solves  the  equilibrium 
boundary-layer  equations  in  a  refined,  near-wall  inner  grid  embedded  in  the  coarser, 
background  LES  grid.  To  model  combustion,  a  simplification  of  the  flamelet-progress 
variable  approach  [FPVA]  of  Pierce  &  Moin  [2004],  with  the  extensions  for  supersonic 
combustion  of  Terrapon  et  al.  [2009,  2010]  and  Pecnik  et  al.  [2012],  was  utilized.  The  FPVA 
model  reduces  the  otherwise  computationally  intractable  complexity  of  a  hydrocarbon  fuel 
chemical  mechanism  [owing  to  the  large  number  of  species  and  reactions  involved]  to  a 
flamelet  look-up  table,  pre-computed  for  a  set  of  flame  boundary  conditions. 
Transport/reaction  equations  are  added  to  the  filtered  Navier-Stokes  equations  in  the 
numerical  solver  for  three  additional  scalar  fields:  the  filtered  mixture  fraction,  Z;  the 
subgrid-scale  variance  of  the  mixture  fraction,  Z"Z";  and  a  filtered  reaction  progress 
variable,  C.  In  the  present  work,  the  progress  variable  is  defined  as  the  sum  of  the  mass 
fractions  of  the  following  combustion  products:  H2O,  CO2,  CO,  and  H2,  and  normalized 
between  0  [non-reacting]  and  1  [fully  reacted].  In  these  simulations,  only  the  fully  reacted 
flamelet  from  the  library  [C  =  1]  is  used.  Note  that  the  flamelet  solution  at  the  boundary 
states  of  pure  oxidizer  and  pure  fuel  coincides  with  the  fully  quenched  [non-reacting] 
flamelet  ((7  =  0).  The  effects  of  turbulence  in  the  subgrid  fluctuations  of  the  mixture  fraction 
and  progress  variable  are  modeled  by  assuming  (3  and  6  probability  density  functions, 
respectively,  in  the  integration  of  combustion  variables  in  ZC-space  leading  to  their  filtered 
counterparts. 

A  volumetric  mixture  of  64%  ethylene  and  36%  methane  is  used  in  the  HIFiRE-2  scramjet 
as  a  two-component  gaseous  surrogate  of  partially-cracked  JP-7,  seeking  to  approximate  its 
ignition  delay,  extinction,  and  flame  strength  characteristics.  The  CEFRC  v.  0.9  mechanism 
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provided  by  the  High  Temperature  Gasdynamics  Laboratory  (HTGDL)  at  Stanford 
University  (CEFRC,  2013)  was  used  to  generate  the  flamelet  library  for  the  simulations 
presented  in  this  brief.  The  mechanism  consists  of  approximately  260  reactions  and  50 
species. 

1.2.2.1  Boundary  conditions 

The  boundary  condition  at  the  inflow  plane  of  the  isolator  matches  the  air  mass  flow  rates 
and  Mach  numbers  given  in  Table  1.  Turbulent  boundary  layers  develop  upstream  of  the 
isolator  along  the  HDCR  nozzle  present  in  the  ground  tests.  We  consider  a  nominal 
boundary  layer  thickness  of  approximately  1  mm  at  the  entrance  of  the  isolator.  We  employ 
a  synthetic  turbulent  inflow  generator  based  on  a  digital  filtering  technique  originally 
proposed  by  Klein  et  al.  [2003],  with  the  improvements  of  Xie  &  Castro  (2008)  and  Touber 
&  Sandham  (2009). 

Characteristic  boundary  conditions  for  velocity,  pressure,  and  temperature  are  used  for  the 
primary  and  secondary  injectors,  matching  the  experimental  fuel  mass  flow  rates  specified 
in  Table  1.  In  the  experiments,  the  fuel  was  heated  to  prevent  liquification  of  the  ethylene 
present  in  the  mixture  (see  Hass  et  al.,  2009).  In  the  simulations,  the  fuel  temperature  is  set 
to  300  K.  The  filtered  mixture  fraction  and  progress  variable  are  set  to  1  (pure  fuel)  and  0 
(non-reacting),  respectively.  The  variance  of  the  subgrid-scale  mixture  fraction  is  set  to 
zero. 

1.2.3  Results 

1.2. 3.1  Time-averaged  pressure  profiles.  Comparison  with  experiments 

The  HDCR  is  instrumented  with  static  pressure  ports  along  the  isolator  and  combustor 
walls.  Fig.  9  shows  a  comparison  of  the  time-averaged  wall  pressure  profiles  obtained  from 
the  simulations  and  the  ground  test  experimental  data  at  the  spanwise  center  plane.  In  the 
simulations,  statistics  are  collected  for  approximately  2  ms  after  initial  transients,  which 
corresponds  to  approximately  4.5  flow-through  times  for  the  flight  Mach  number  of  6.5, 
based  on  the  centerline  velocity  at  the  entrance  to  the  isolator. 

For  the  dual-mode  operation  at  Mf  =  6.5,  the  simulations  are  started  from  quiescent  flow  in 
the  non-reacting  regime  (without  fuel  injection).  The  hollow  symbols  and  dash-dotted  line 
in  Fig.  9(a)  correspond  to  the  experiments  and  (coarse-mesh)  simulation  results  for  that 
non-reacting  (tare)  case.  The  shape  of  the  pressure  profile  is  well  captured  by  the 
simulation.  The  pressure  recovery  along  the  cavity  ramp  appears  higher  in  the  simulations, 
which  might  be  a  result  of  differences  in  the  reattachment  location  with  respect  to  the 
experiments  associated  with  the  wall  model. 

Results  with  the  flamelet-based  fast-chemistry  combustion  model  are  shown  for  the  three 
grid  resolutions  under  consideration  (in  dotted,  dashed,  and  solid  lines,  for  increasing 
resolution).  Experimental  pressure  levels  inside  the  cavity  and  in  the  last  two  thirds  of  the 
secondary  combustor  are  recovered  in  the  simulations.  The  pressure  rise  upstream  of  the 
primary  injector  is  also  captured  although  it  is  found  downstream  of  its  corresponding 
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experimental  location,  which  penetrates  farther  into  the  isolator.  A  similar  outcome  was 
found  by  Bynum  &  Baurle  (2011)  who  performed  an  uncertainty  quantification  study  of  the 
variability  of  the  HIFiRE-2  ground  test  for  the  dual-mode  regime  (at  a  lower  Mf  =  5.84).  It  is 
worth  noting  that  these  simulations  do  not  currently  account  for  any  possible  incoming 
nitric  oxide  present  in  the  incoming  air  resulting  from  the  arc-heating  (Hass  etal.,  2009), 
which  might  be  responsible  for  enhancing  combustion  (Pellett  etal.,  2009),  and  could 
contribute  to  the  combustion-induced  pressure  rise  seen  to  occur  experimentally  farther 
upstream  in  the  isolator. 


Fig.  9.  Time-averaged  wall-pressure  profiles  at  the  center  plane,  (a)  Dual-mode  regime  at 
Mf  =  6.5:  symbols  correspond  to  the  experiment  (hollow  for  non-reacting  and  solid  for 
reacting),  and  lines  correspond  to  simulation  results:  dash-dotted  for  non-reacting;  dotted, 
dashed  and  solid  for  reacting  with  flamelet-based  fast-chemistry  combustion  model  on 
coarse-,  medium-,  and  fine-resolution  meshes,  respectively,  (b)  Transition  from  dual-  to 
scram-mode:  symbols,  experimental  data  at  Mf  =  8.0;  lines,  WMLES  results  with  flamelet- 
based  fast-chemistry  combustion  model  on  medium-resolution  mesh,  with  each  line 
corresponding  to  the  simulation  data  averaged  over  consecutive  time  intervals  of  0.3  ms, 
evolving  from  the  sudden  change  of  conditions  from  dual-  to  scram-mode.  For  reference,  a 
side  view  of  the  geometry  of  the  scramjet  is  shown  on  top  of  each  plot,  along  with  vertical 
dotted  lines  indicating  the  primary  and  secondary  stages  of  fuel  injection. 

To  simulate  the  dual  to  scramjet-mode  transition,  from  the  solution  of  the  dual-mode 
simulations  at  Mf  =  6.5,  a  sudden  change  in  the  boundary  conditions  at  the  entrance  to  the 
isolator  and  the  fuel  injectors  is  imposed  to  match  the  scramjet  operation  at  Mf  =  8  in  the 
ground  tests.  The  flow  then  evolves  through  a  transient  reflected  in  Fig.  9(b),  which  shows 
center-line  wall-pressure  profiles  averaged  over  consecutive  time  intervals  of  0.3  ms.  It  is 
observed  that  within  the  first  2  ms,  the  pressure  levels  throughout  the  scramjet  decrease  to 
the  experimental  values. 

Fig.  10(a-g)  shows  contours  of  instantaneous  thermodynamic  variables,  species 
concentration  and  velocity  fields  extracted  on  vertical  and  horizontal  slices  of  the  flow 
domain.  The  magnitude  of  the  density  gradient  on  the  vertical  plane  shows  the  oblique 
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shock  generated  at  the  end  of  the  isolator,  whereas  the  temperature  field  confirms  that 
mixing  and  combustion  occur  in  a  recirculation  region  upstream  of  the  primary  injector. 
The  concentration  of  species  for  different  combustion  products  in  Fig.  10(c,  d  and  h)  also 
confirms  the  combustion-induced  origin  of  the  pressure  increase  upstream  of  the  primary 
injection.  Peak  values  of  OH  are  seen  near  the  wall  in  the  cavity  ramp,  where  the  wall- 
pressure  levels  are  higher  (see  Fig.  9a]  and  the  mixture  fraction  is  closer  to  its 
stoichiometric  value  (enveloped  by  the  black  lines  in  Fig.  lOe).  In  contrast,  H2O  shows  a 
more  homogeneous  concentration  inside  the  first  part  of  the  cavity  (i.e.,  the  part  with  the 
nearly  constant  cross  section),  decreasing  along  the  cavity  ramp;  CO  concentration  (shown 
in  Fig.  lOh  for  a  plane  parallel  to  the  combustor  wall  separated  2  mm  from  it)  peaks  first 
along  the  initial  part  of  the  combustor,  upstream  of  the  cavity.  The  mixture  fraction 
(Fig.  lOe)  appears  relatively  more  uniform  inside  the  cavity  (where  the  maximum 
temperatures  are  reached)  than  downstream  of  the  secondary  injection,  where  large-scale 
unsteadiness  is  present  and  results  in  more  complex  mixing  and  combustion  patterns, 
which  are  also  observed  in  the  contours  of  the  density-gradient  magnitude  and  the 
concentration  of  combustion  products. 


Fig.  10.  From  top  to  bottom:  contours  of  instantaneous  temperature  (215-3280  K),  density 
gradient  (15-500  kg/m4,  logscale),  OH  concentration  (0-0.01),  H2O  concentration  (0-0.09), 
mixture  fraction  (0-1),  pressure  (77-773  KPa),  streamwise  velocity  (-600-1615  m/s),  and 
CO  concentration  (0-0.3)  for  the  reacting-flow  simulation  in  the  dual-mode  operating 
regime  at  Mf  =  6.5.  The  bottom  plot  (h)  contains  a  slice  of  the  computational  domain 
parallel  to  the  wall  of  the  combustor  at  a  distance  2  mm  normal  to  the  wall.  All  other  plots 
(a-g)  are  vertical  slices  at  z  =  12.7  mm  passing  through  a  set  of  injectors.  For  clarity  most  of 
the  isolator  has  been  left  out  of  the  visualization. 
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Fig.  11  shows  a  comparison  between  simulation  results  in  the  dual-  and  scram-mode 
operating  regimes  at  Mf  =  6.5  and  8.0,  respectively.  The  instantaneous  Mach  number 
contours  (Fig.  19a]  confirm  that  in  the  scram-mode,  the  flow  remains  fully  supersonic 
along  the  core  of  the  engine,  with  the  subsonic  flow  regions  mostly  confined  to  the  cavity 
and  some  patches  in  the  wake  of  the  jet  downstream  of  the  secondary  injector.  In  the  dual¬ 
mode,  the  shock-train  inside  the  combustor  extends  upstream  near  the  isolator,  with 
subsonic  flow  dominating  over  a  larger  portion  of  the  combustor  and  extending  to  the  core 
of  the  engine. 


Fig.  11:  Contours  of  instantaneous  Mach  number  (a]  and  CO  concentration  (b)  on  the 
spanwise-normal  plane  located  at  z  =  12.7  mm.  Each  plot  includes  the  results  from  the 
scram-mode  (Mf  =  8]  and  dual-mode  (Mf  =  6.5)  simulations  plotted  on  the  half  top  and 
bottom  parts,  respectively,  using  the  engine  symmetry,  for  comparison.  Colorbar  for  the 
Mach  number  from  0  to  3.6  (with  the  sonic  line  in  white).  Colorbar  for  Y  co  from  0  to  0.3. 

The  recirculation  region  upstream  of  the  primary  injector  for  the  dual-mode  is  not  present 
in  the  scram-mode.  The  contours  of  CO  concentration  in  Fig.  11(b)  suggest  that  the 
combustion  of  the  fuel  injected  in  the  primary  stage  occurs  rather  differently  between  the 
two  modes:  in  the  dual-mode,  the  highest  CO  concentration  levels  occur  shortly 
downstream  of  the  primary  injector,  decreasing  inside  the  cavity;  however,  in  the  scram- 
mode  simulation,  only  a  thin  layer  of  CO  is  visible  immediately  downstream  of  the  primary 
injector,  with  most  of  the  CO  concentrating  inside  the  cavity  flameholder.  Downstream  of 
the  secondary  injector,  the  CO  concentration  levels  and  patterns  appear  similar  between 
the  two  modes. 


Part  2 

Large  Eddy  Simulation  of  Crackle  Noise  in  Strongly 

Heated  Supersonic  Jets 


2.1  Background 

Crackle,  first  investigated  by  Ffowcs  Williams  et  al.  (1975),  is  a  component  of  supersonic  jet 
noise,  which  is  particularly  irritating  when  it  occurs.  Crackle  is  characterized  by 
intermittent  positive  pressure  "spikes"  caused  by  N-shaped  waves  arriving  at  observer 
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locations.  Crackle  radiation  is  directed  downstream  at  an  angle  associated  with  the  peak  jet 
noise,  and  may  account  for  as  much  as  30%  of  the  overall  sound  pressure  level  in  this 
direction  (Krothapalli  et  al.,  2000,  2003).  Therefore,  the  elimination  of  crackle  has  the 
potential  to  reduce  peak  jet  noise  levels  by  3-5dB,  which  would  meet  near-term  jet  noise 
reduction  goals  (Martens  &  Spyropoulos  2010). 

The  mechanism  by  which  N-shaped  waves  are  generated  is  not  yet  fully  understood.  In 
particular,  there  is  debate  about  whether  these  waves  are  a  result  of  nonlinear  acoustic 
propagation  or  whether  they  are  generated  inside  the  supersonic  jet  itself.  Recent  direct 
numerical  simulations  of  temporal  supersonic  shear  layers  show  wave  agglomeration,  a 
nonlinear  effect  where  strong  acoustic  waves  travel  slightly  faster  than  weaker  waves 
(Andersen  &  Freund  2012).  The  strong  waves  accumulate  the  energy  of  the  weaker  waves 
as  they  overtake  them.  While  simulations  of  a  plane  shear  layer  may  model  the  flow  near 
the  nozzle  lip,  such  an  approach  necessarily  neglects  large-scale  effects  associated  with  the 
cylindrical  nature  of  the  jet.  Such  effects  include  instability  waves  for  which  azimuthal 
wavenumbers  m  =  0  (axisymmetric)  and  m  =  1  (helical)  are  the  most  prevalent. 
Furthermore,  simulations  of  a  temporal  shear  layer  with  streamwise  periodic  boundary 
conditions  neglects  effects  of  the  spatial  development  of  the  shear  layers.  The  instability  of 
a  temporal  simulation  gradually  feeds  energy  to  increasingly  large  scales  over  time, 
whereas  such  growth  is  convected  downstream  in  a  spatially  developing  layer.  In  a 
temporally  evolving  shear  layer,  the  agglomeration  of  acoustic  waves  over  time  could  be  a 
consequence  of  this  nonphysical  effect. 

On  the  other  hand,  recent  experiments  have  provided  some  evidence  that  crackle  may  be  a 
property  of  the  aerodynamics  within  the  jet  itself.  For  example,  it  was  found  that  the 
addition  of  chevrons  to  the  lip  of  a  military-style  nozzle  significantly  reduces  crackle  noise 
emissions.  Likewise,  twin-notched  nozzles  and  nozzles  with  efficient  convoluted  silencers 
are  also  found  to  reduce  crackle  (Ffowcs  Williams  et  al.  1975).  The  sensitivity  of  crackle  to 
nozzle  lip  modifications,  which  affect  the  issuing  jet  plume,  suggest  that  the  N-shaped 
waves  perceived  as  crackle  are  generated  in  the  jet  itself.  The  aerodynamic  mechanism  by 
which  supersonic  jets  create  these  waves  remains  a  mystery.  Our  approach,  discussed 
below,  is  to  complement  experimental  measurements  by  simulating  an  entire  spatially 
developing  heated  supersonic  jet  issuing  from  a  realistic  military-style  nozzle  using  high- 
fidelity  Large  Eddy  Simulation  (LES)  (Nichols  et  al.  2013).  To  meet  this  challenge  of 
simulating  a  broad  range  of  length  scales,  we  apply  LES  on  unstructured  meshes  to  allow 
grid  points  to  be  efficiently  clustered  where  needed.  The  unstructured  solver  CharLES, 
developed  at  Cascade  Technologies,  Inc.,  efficiently  simulates  turbulence  on  unstructured 
meshes  by  using  widened  numerical  stencils  to  explicitly  control  dissipation.  This  method 
has  been  extensively  tested  in  the  context  of  the  prediction  of  the  aeroacoustics  of  jets  from 
complex  geometries  (Nichols  et  al.  2011).  The  seamless  treatment  of  complex  geometry  is 
particularly  important  when  nozzle  lip  modifications  are  considered  for  crackle 
suppression. 
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2.2  Approach 


As  a  numerical  test  case,  we  investigate  crackle  produced  by  the  military-style  nozzle 
shown  in  Fig.  1(a)  (General  Electric  nozzle  L03 116-410).  This  nozzle  incorporates  a  conical 
contraction  connected  to  a  faceted  straight-ramp  diffuser  at  a  relatively  sharp  throat.  In  the 
diffuser  section,  the  twelve  facets  are  created  by  what  would  be  twelve  seals  in  an  actual 
variable  area  engine,  although  in  the  experimental  model,  the  facets  are  fixed.  The  area 
ratio  of  the  nozzle  exit  to  the  nozzle  throat  is  1.295,  so  that  the  design  Mach  number  is 
Md  =  1.65.  The  simulation  discussed  in  the  following  text  was  driven  by  a  nozzle  pressure 
ratio  of  NPR  =  P0/Pm  =  4.0  and  a  nozzle  temperature  ratio  of  NTR  =  T0/Tw  =  3.65,  where 
the  subscript  0  refers  to  stagnation  properties  inside  the  nozzle  and  the  subscript  oo  refers 
to  static  properties  in  the  ambient  fluid.  This  operating  point  corresponds  to  a  heated 
supersonic  jet  with  a  fully-  expanded  Mach  number  of  M;-  =  Uj/cj  =  1.56  and  an  acoustic 
Mach  number  of  Ma  =  Uj/cm  =  2.44,  where  Uj  is  the  jet  velocity,  and  Cj  is  the  speed  of 
sound  within  the  jet  plume,  which  is  greater  than  the  speed  of  sound  cm  in  the  cooler 
ambient  fluid.  Because  the  fully  expanded  jet  Mach  number  is  less  than  the  design  Mach 
number,  the  flow  is  over-expanded  as  it  leaves  the  nozzle  and  produces  a  train  of  shock 
cells  downstream.  The  over-expanded  condition  was  chosen  since  the  normal  condition  at 
take-off  also  involves  over-expansion.  To  further  simulate  the  conditions  at  initial  takeoff, 
no  coflow  was  used  in  the  simulations. 


(a) 


Fig.  12.  (a)  A  military-style  nozzle  designed  by  GE  with  a  faceted  straight-ramp  diffuser 
used  for  the  crackle  simulations,  (b)  An  axial  cross  section  showing  the  mesh  resolution  in 
the  near-nozzle  region.  A  zone  extending  8  jet  diameters  downstream  of  the  jet  exit  was 

adaptively  refined. 


Taking  advantage  of  the  unstructured  mesh  capability  provided  by  CharLES,  a  body-fitted 
mesh  was  generated  precisely  conforming  to  all  of  the  geometric  details  of  the  nozzle, 
including  the  facets,  the  small  radius  at  the  nozzle  throat,  the  nozzle  lip,  and  the  shroud 
housing  the  nozzle  assembly.  A  section  of  the  nozzle  facility  adapter  was  also  incorporated 
into  the  computational  domain  in  order  to  allow  the  flow  injected  at  the  upstream 
boundary  to  naturally  develop  before  encountering  the  contraction.  To  capture  the  jet 
plume,  the  computational  domain  was  extended  a  distance  of  27.5  nozzle  diameters 
downstream  of  the  nozzle  exit.  A  mesh  ideally  suited  to  capture  acoustic  production  (see 
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Nichols  et  al.  2011]  was  generated  by  first  constructing  a  coarse  (~6xl06  control  volumes) 
body-fitted  mesh  and  then  applying  an  adaptive  refinement  procedure  so  that  the 
turbulence  containing  regions  in  the  jet  plume  are  resolved  by  a  nearly  isotropic  mesh.  An 
additional  level  of  refinement  was  used  for  the  near  nozzle  region  to  capture  the  evolution 
of  the  jet  shear  layers  as  they  emerge  from  the  nozzle.  Figure  1(b)  provides  a  visualization 
of  the  initial  portion  of  the  resulting  mesh  which  contained  ~331xl06  hexahedral  control 
volumes.  Characteristic  boundary  conditions  were  implemented  at  all  of  the  boundaries, 
accompanied  by  numerical  sponge  layers  designed  to  minimize  unphysical  acoustic 
reflection.  In  particular,  the  outlet  sponge  was  carefully  designed,  following  the  guidelines 
put  forth  by  Mani  (2012). 

2.3  Results 

Figure  13  shows  an  axial  cross-section  through  a  snapshot  of  the  temperature  and  pressure 
fields  taken  from  the  LES  of  the  heated  supersonic  jet.  In  the  figure,  yellow  scale  contours  of 
temperature  are  shown  in  the  interior  of  the  jet,  and  show  the  turbulent  eddies  within  the 
jet  plume.  Also  visible  in  the  temperature  field  is  a  train  of  shock  cells,  initiated  for  this 
straight-ramp  diffuser  at  both  the  sharp  nozzle  throat  and  the  nozzle  lip,  leading  to  a 
double-diamond  pattern  downstream.  While  the  jet  is  over-  expanded,  the  pressure 
mismatch  is  not  too  severe,  thus  the  shocks  associated  with  these  cells  are  weak.  Exterior 
to  the  jet,  blue  scale  pressure  contours  are  used  to  visualize  the  corresponding  near-field 
acoustics  generated  by  the  jet.  Both  the  temperature  and  pressure  have  been 
nondimensionalized  by  the  ambient  conditions.  The  figure  shows  that  the  near-field 
acoustics  are  dominated  by  a  downstream  component.  Upstream-propagating  components 
corresponding  to  broadband  shock-associated  noise  can  also  be  observed.  These  are 
relatively  weak,  however,  indicating  that  the  jet  is  not  too  far  from  ideal  expansion. 

The  geometry  shown  in  Fig.  12(a)  exactly  matches  that  of  an  experimental  nozzle  used  for 
acoustic  testing  (Martens  et  al.  2011).  To  validate  our  numerical  methodology,  we  use  the 
LES  data  to  compute  the  farfield  noise  spectrum  at  the  same  location  as  measured  from  the 
experiment.  For  this  purpose,  the  acoustic  data  is  recorded  along  a  surface  surrounding  the 
turbulent  jet  plume  and  then  propagated  to  the  farfield  by  applying  the  Ffowcs  Williams- 
Hawkings  (FWH)  equation.  Unfortunately,  because  the  FWH  equation  is  based  on  an 
assumption  of  linearity,  we  do  not  expect  this  method  to  be  able  to  exactly  reproduce  the 
crackle  noise  in  the  farfield.  The  FWH  equation,  by  definition,  neglects  the  nonlinear 
acoustic  propagation  effects  that  may  alter  the  N-shaped  waveforms  over  a  distance, 
although  for  the  amplitudes  considered,  we  expect  this  influence  to  be  relatively  small. 
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Fig.  13.  An  instantaneous  snapshot  of  the  heated  supersonic  jet  visualized  by  contours  of 
temperature  [interior  of  the  jet,  yellow  scale)  and  contours  of  pressure  (exterior  of  the  jet, 
blue  scale)  on  an  axial  cross  section.  The  nozzle  cross-section  appears  in  green  at  left. 

Figure  14  compares  the  farfield  acoustic  spectrum  measured  from  the  experiment  (black 
curve)  to  that  computed  from  the  simulation  data  using  the  FWH  equation  (red  curve).  The 
measurement  location  was  at  an  angle  of  140  degrees  with  respect  to  the  upstream  jet  axis 
and  at  a  distance  of  72 D  from  the  center  of  the  nozzle  exit.  This  corresponds  to  the 
direction  of  the  peak  jet  noise.  In  the  figure,  the  blue  curve  represents  the  narrowband 
spectrum  computed  from  the  simulation  to  which  frequency  binning  was  applied  to 
compute  the  third  octave  spectrum  (red)  at  the  same  frequencies  reported  by  the 
experiment.  Although  there  are  some  discrepancies  between  the  experiment  and 
simulation  data,  the  overall  shape  and  levels  of  the  curve  agree  reasonably  well  over  a 
broad  range  of  frequencies,  giving  us  confidence  in  the  accuracy  of  our  simulation. 


Fig.  14.  Acoustics  at  an  angle  of  140  deg  to  the  jet  upstream  axis  (peak  jet  noise  direction). 
A  comparison  of  the  experimental  and  simulated  farfield. 
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In  the  near  field  of  the  jet,  Mach  wave  radiation  at  an  angle  of  approximately  140  degrees  is 
clearly  evident  in  Fig.  13.  In  addition,  interspersed  throughout  the  Mach  wave  field,  we 
observe  thin  structures  of  intense  pressure  in  the  near  field  of  the  jet.  In  these  structures, 
the  pressure  perturbation  is  about  10%  of  the  ambient  pressure,  corresponding  to  a  strong 
acoustic  wave  or  a  weak  shock.  We  also  note  that  these  shocklet  structures  consist  of  a 
sharp  front  immediately  followed  by  a  more  gentle  decay  in  pressure.  As  a  shocklet 
propagates  outwards,  away  from  the  jet,  it  creates  a  sudden  compression  followed  by  a 
more  gradual  expansion  at  any  fixed  observer  location,  which  can  then  be  interpreted  as  a 
crackle  event. 

To  emulate  the  near-field  acoustic  array  used  in  previous  scale-model  crackle  experiments 
(Martens  et  al.  2011],  we  probe  the  simulation  at  the  locations  indicated  by  the  circles  in 
Fig.  13.  These  points  correspond  to  sixteen  unique  axial  and  radial  positions  arranged  into 
inner  and  outer  arrays  of  eight  points  each.  For  each  of  these  probes,  data  was  taken  at  48 
equidistant  locations  along  a  ring  in  the  azimuthal  direction.  In  the  simulation,  we  angle  the 
near-field  probe  array  with  the  jet  spreading,  instead  of  collecting  data  at  a  uniform  radial 
distance  away  from  the  jet,  as  was  done  in  the  laboratory  experiment.  This  was  done 
because  the  quality  of  the  mesh  rapidly  degrades  as  we  move  away  from  the  turbulence 
containing  region;  by  remaining  a  constant  distance  from  the  turbulence,  we  ensure  a 
constant  mesh  quality  across  all  of  the  probes.  The  simulation  environment  allows  us  to 
determine,  with  great  precision,  exactly  where  the  turbulent  region  ends  so  that  we  may 
accurately  place  a  probe  very  close  (much  closer  that  would  be  safe  in  experiments)  to  the 
turbulence  in  order  to  capture  the  acoustic  field  as  it  first  emerges  from  the  turbulence. 
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Fig.  15.  Pressure  signal  at  inner  probe  location  3,  with  skewness  level  0.4396.  The  dashed 
lines  are  spaced  at  standard  deviations  away  from  the  mean  (red  line).  Crackle  is 
characterized  by  large  positive-pressure  excursions  and  N-shaped  waves. 

Figure  15  shows  a  time  history  of  the  local  pressure  recorded  at  the  third  inner  probe 
location.  The  simulation  was  run  for  approximately  450  nondimensional  time  units  based 
on  the  characteristic  time  D/uj.  The  horizontal  solid  red  line  indicates  the  signal  mean  and 
the  dashed  lines  are  spaced  at  one  standard  deviation.  The  signal  shown  in  Fig.  15  is 
skewed  towards  the  positive  side  and  makes  large  intermittent  positive-pressure 
excursions,  reaching  five  standard  deviations  away  from  the  mean  in  some  cases.  On  the 
negative  pressure  side,  the  signals  rarely  venture  even  three  standard  deviations  away 
from  the  mean.  The  skewness  s  of  the  pressure  signals  is  defined  as 
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E[(.P(t)  -  Li)3] 

s  = - 3 - 

as 

where  jU  is  the  signal  mean  and  a  is  the  standard  deviation.  The  skewness  of  the  signal 
shown  in  Fig.  15  is  0.4396,  which  meets  the  criterion  that  s  >  0.4  indicates  crackle,  as  put 
forth  by  Ffowcs  Williams  et  al.  [1975).  Furthermore,  to  investigate  the  spatial  dependence 
of  the  skewness,  we  compute  the  skewness  of  the  signals  at  each  of  the  48  azimuthal  probe 
locations  and  then  take  an  average.  Fig.  16  shows  the  results,  along  with  error  bars 
indicating  the  standard  deviation  of  the  levels  over  the  48  probe  locations.  In  the  figure,  the 
solid  line  corresponds  to  the  inner  probe  array  shown  in  Fig.  13,  while  the  dashed  line 
corresponds  to  the  outer  probe  array.  The  skewness  is  observed  to  have  a  maximum  close 
to  the  nozzle  and  another  maximum  farther  downstream  [near  x/D  =  10). 


Fig.  16.  Skewness  levels  versus  axial  position  along  the  inner  [solid  curve)  and  outer 
[dashed  curve)  probe  arrays.  At  each  station,  the  skewness  was  averaged  with  respect  to 
the  azimuthal  direction.  The  error  bars  indicate  standard  deviation  of  the  azimuthal  data. 

Figure  17[a)  visualizes  the  skewness  along  a  portion  of  the  FWH  surface  used  to  compute 
the  spectra  shown  in  Fig.  14.  The  portion  of  the  FWH  surface  shown  has  the  shape  of  a 
truncated  cone  and  extends  from  the  plane  of  the  nozzle  exit  to  approximately  eight 
diameters  downstream.  This  corresponds  to  the  portion  of  the  FWH  surface  that  passes 
through  the  near-  nozzle  refinement  region  shown  in  Fig.  12[b).  On  this  surface,  contours 
of  pressure  are  shown.  Figure  17[b)  shows  the  same  surface  as  in  Fig.  17[a),  but 
unwrapped,  so  the  horizontal  axis  is  the  axial  distance,  while  the  vertical  axis  gives  the 
azimuthal  angle.  Figure  17[b)  shows  that  the  waveforms  associated  with  the  highest  levels 
of  skewness  tend  to  be  correlated  over  60  degrees  of  the  azimuthal  extent,  although  much 
smaller  patches  of  high  pressure  also  appear.  Note  also  that  each  patch  of  high  skewness 
forms  a  sharp  gradient  along  its  downstream  edge,  indicating  an  N-shaped  waveform. 

To  further  understand  the  source  of  the  crackle,  we  consider  the  unsteady  flow  evolution 
leading  up  to  a  typical  crackle  event.  For  this  purpose,  we  consider  the  outer  probe 
location  2  where  the  skewness  is  measure  to  be  0.4245.  Figures  18[a)  provides  a  close  up 
view  of  the  pressure  signals  at  times  near  to  the  crackle  event.  Note  that  the  crackle  event 
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corresponds  to  a  sharp  compression  followed  by  a  gradual  relaxation.  Such  N-shaped 
waveforms  are  also  evident  in  the  pressure  signals  measured  along  the  inner  array.  The 
fact  that  such  waveforms  are  fully  formed  this  close  to  the  jet  strongly  supports  the  notion 
that  crackle  waves  are  formed  directly  in  the  jet  rather  than  as  a  consequence  of  nonlinear 
propagation. 

(a)  (b) 
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Fig.  17.  (a)  Pressure  contours  along  a  portion  of  the  conical  FWH  surface,  (b)  The 
“unwrapped"  surface  showing  azimuthal  dependences  of  skewed  pressure  waves  along  the 

entire  circumference. 
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Fig.  18.  (a]  Pressure  signal  at  times  near  a  typical  crackle  event  (marked  by  arrow),  (b)- 
(d)  Time  sequence  of  the  flow  field  evolution  leading  up  to  this  event.  The  jet  interior  and 
shears  layers  are  visualized  by  contours  of  temperature  (grayscale),  while  pressure 
contours  are  shown  exterior  to  the  jet  (color).  The  solid  magenta  curve  represents  the 
sonic  line,  which  separates  fluid  moving  supersonically  with  respect  to  the  ambient  from 
fluid  moving  subsonically.  The  circle  indicates  the  probe  location  and  the  arrows  indicate 
the  wave  that  causes  the  crackle  event  at  the  final  time. 

Figures  18(b)-18(d)  visualize  the  entire  flow  field  at  the  instants  leading  up  to  the  crackle 
event.  In  these  figures,  the  probe  location  (outer  probe  2)  is  indicated  by  the  green  circle. 
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Similar  to  Fig.  13,  temperature  contours  are  shown  for  the  jet  interior  and  shear  layers, 
while  pressure  contours  are  shown  exterior  to  the  jet.  The  pressure  contours  use  the  same 
color  scale  as  Fig.  17,  to  highlight  pressures  more  than  three  standard  deviations  above  the 
mean.  The  solid  magenta  curve  indicates  the  contour  where  the  acoustic  Mach  number  is 
equal  to  unity.  Inside  this  contour,  the  fluid  is  moving  supersonically  with  respect  to  the 
ambient  fluid,  while  outside  the  fluid  is  moving  subsonically. 

In  Fig.  18(b),  we  observe  that  the  creation  of  the  crackle  wave  coincides  with  a  large  bulge 
of  the  sonic  line  just  upstream.  From  the  temperature  contours,  the  foot  of  the  emerging 
crackle  wave  appears  to  be  anchored  in  the  region  between  two  eddies.  In  this  region,  a 
tongue  of  cold  ambient  fluid  (darker  gray)  is  being  entrained  into  the  hot  jet  core  by  the 
action  of  these  two  vortices.  The  entrainment  of  cold  fluid  causes  the  sonic  line  to  slightly 
indent  in  this  region.  Also,  because  this  interface  is  being  stretched,  a  sharp  front  forms 
between  the  hot  and  cold  fluids.  This  pattern  is  also  evident  at  the  foot  of  the  companion 
crackle  wave  just  downstream  (which  shows  up  later  at  probe  2  as  the  pressure  spike  just 
before  the  crackle  event).  In  Fig.  18(c),  we  see  that  the  indentation  of  the  sonic  line  has 
propagated  downstream,  along  with  the  foot  of  the  crackle  wave.  If  the  propagation  of 
deflections  of  the  sonic  line  propagates  supersonically  with  respect  to  the  ambient  fluid, 
then  we  may  apply  the  supersonic  wavy  wall  analogy  (Tam  1995)  leading  to  Mach  waves. 


Fig.  19.  Contours  of  the  shock  sensor  C  on  an  axial  cross  section  through  the  heated  jet, 
highlighting  shocklets  in  the  jet  near  field.  Dashed  lines  delimit  regions  of  near-nozzle 
adaptive  refinement.  Circles  indicate  locations  where  the  shocklet  foot  is  embedded  inside 

the  turbulent  shear  layer. 

Finally,  in  order  to  examine  the  shocklet  structures  in  the  flow,  we  apply  a  diagnostic 
inspired  by  the  shock  sensor  of  Ducros  et  al.  (1999),  modified  by  Bhagatwala  and  Lele 
(2009,2011)  to  highlight  regions  of  negative  dilatation  with 
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where  fi  is  the  vorticity  magnitude,  D  is  the  grid  spacing,  and  c  is  the  local  speed  of  sound. 
Figure  19  shows  contours  of  this  quantity  at  a  single  time  instant.  The  dashed  lines  delimit 
the  near-nozzle  region  of  mesh  refinement  as  shown  in  Fig.  13.  In  the  core  of  the  jet,  the 
shocks  from  the  sharp  nozzle  throat,  the  nozzle  lip,  and  the  resulting  shock  cells 
downstream  are  clearly  visible  using  the  shock  sensor  C.  The  turbulent  shear  layers  of  the 
jet  are  highly  vortical,  so  that  the  term  H2  in  the  denominator  causes  C  to  be  small.  In  Fig.  8, 
the  turbulent  shear  layers  appear  as  shadows  immediately  above  and  below  the  jet  core. 
Exterior  to  the  turbulence,  the  shock  sensor  highlights  many  thin  regions  of  strong 
compression  emanating  from  the  jet.  These  regions  of  negative  dilatation  correspond  to  the 
sharp  fronts  of  the  N-shaped  waves,  or  shocklets,  emitted  from  the  jet.  While  the  shocklets 
eventually  detach  from  the  jet  and  propagate  away  to  the  far  field,  Fig.  19  shows  several 
locations  (highlighted  by  circles)  where  the  foot  of  an  emerging  shocklet  is  embedded 
directly  inside  the  turbulent  shear  layer.  At  these  locations,  we  note  that  the  vorticity 
magnitude  must  not  be  too  large;  otherwise,  C  would  be  attenuated.  The  embedded 
shocklets  appear,  instead,  to  be  formed  in  the  interstices  between  two  regions  of  higher 
vorticity:  the  braid  region  between  two  corotating  large-scale  coherent  fluid  motions,  see 
for  example  Nichols  et  al.  (2014). 

2.4  Crackle  noise  conclusions 

An  unstructured  LES  of  a  turbulent  supersonic  jet  issuing  from  a  faceted  military-style 
nozzle  was  performed  at  operating  conditions  NPR  =  4.0  and  NTR  =  3.65,  which  produced 
a  crackling  jet.  Skewness  levels  were  measured  from  the  direct  near  acoustic  field  at 
various  axial  and  radial  positions.  It  was  found  that  the  pressure  perturbations  were  highly 
skewed  very  close  to  the  jet.  Furthermore,  we  observed  crackle  waves  emerging  directly 
from  the  jet  turbulence,  with  characteristic  N-shaped  pressure  waveforms.  Therefore, 
while  nonlinear  propagation  effects  may  eventually  steepen  waves  yet  further,  it  is  not  a 
necessary  component,  since  steep  crackle  waves  are  produced  directly  at  the  jet.  Crackle 
waves  are  produced  at  indentations  of  the  sonic  line  where  cool  ambient  fluid  is  entrained 
into  the  jet  core. 

Part  3 

Linear  stability  analysis  of  the  onset  dynamics  of 

scram  jet  unstart 

3.1  Motivation  and  objectives 

Recently,  interest  has  increased  in  using  scramjet  engines  as  a  means  of  long-range  high¬ 
speed  flights  and  economical  access  to  outer  space.  One  of  the  most  perilous  causes  of 
scramjet  malfunctions  is  the  unstart  event  that  is  initiated  by  excessive  heat  release  from 
combustion.  When  unstart  occurs,  a  strong  moving-shock  structure  is  first  formed  in  the 
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engine,  and  the  shock  structure  propagates  upstream  and  finally  spills  out  of  the  engine 
inlet.  The  unstart  event  is  detrimental  to  the  engine  because  [1]  the  moving-shock 
structure  imposes  high  pressure  and  thermal  loads  on  the  inner  walls  of  the  engine  during 
the  unstart  process,  and  the  walls  can  be  ruptured  due  to  the  loads;  and  [2]  the  airflow  into 
the  engine  is  greatly  diminished  when  the  shock  structure  is  disgorged  by  the  engine, 
leading  to  loss-of-thrust  and  engine  stall.  Because  the  probability  of  unstart  grows  with 
increasing  heat  release  from  combustion,  the  danger  of  unstart  is  an  important  limiting 
factor  in  the  performance  of  scramjets. 

Therefore,  the  onset  mechanisms  of  the  unstart  event  need  to  be  understood  to  prevent  or 
delay  the  unstart  process.  However,  the  detailed  dynamics  has  not  been  fully  understood 
yet,  even  though  many  studies  have  examined  unstart  onset  mechanisms.  For  instance, 
Korkegi  [1975]  suggested  empirical  correlation  functions  for  estimating  the  critical 
pressure  rise  above  which  unstart  occurs,  based  on  the  assumption  that  shock-induced 
flow  separation  of  turbulent  boundary  layers  on  the  engine  walls  causes  the  unstart 
process.  In  ground  tests  of  the  HyShot  II  scramjet  model  (Frost  et  al.,  2009],  the  critical 
pressure  rise  in  the  model  agreed  with  the  Korkegi  limit,  and  therefore  the  authors 
presumed  that  unstart  was  initiated  by  flow  separation  of  the  boundary  layers.  In  a  later 
ground  experiment  of  HyShot  II,  however,  Laurence  etal.  (2013)  could  not  find  large-scale 
boundary-layer  separations,  and  they  concluded  that  flow  separation  was  not  the  main 
cause  of  unstart.  Instead,  they  proposed  thermal  choking  as  the  responsible  unstart 
mechanism.  However,  further  conclusions  regarding  the  onset  mechanism  could  not  be 
drawn  because  of  the  limited  diagnostics  in  the  experiments. 

The  primary  objective  in  this  study  is  to  find  the  onset  dynamics  of  the  unstart  event  based 
on  linearized  system.  Section  2  describes  the  linearized  system  dynamics  that  will  be 
discussed  throughout  this  study.  In  section  3,  the  detailed  methodology  and  the  scramjet 
configuration  are  presented.  The  linearized  dynamics  at  the  unstart  onset  point  is  then 
discussed  in  section  4.  Finally,  section  5  summarized  the  findings  and  make  suggestions  for 
future  work. 


3.2  Linearized  dynamics  of  a  scramjet  system 


In  this  study,  we  consider  the  dynamic  system  of  the  following  form: 


dw 

dt 


F(w;  (fi),  w—[p.f)U,pE,  scalars 


where  c|)  is  an  adjustable  input  parameter  to  the  system,  w  is  a  solution  vector  describing 
the  system  state,  and  F  is  a  non-linear  map  describing  the  temporal  evolution  of  the  system 
state.  In  this  study,  we  consider  supersonic  flows  in  a  model  scramjet  with  a  heat  source. 
Therefore,  F  is  the  compressible  Navier-Stokes  equations  with  a  turbulent  model,  and  4>  is 
the  heat  release  rate  of  the  source,  w  is  the  flow  variable  vector  whose  elements  are  density 
(p),  the  momentum  vector  (pu  =  [pui,pu2,pu3]T],  total  energy  (pE),  and  conservative  scalars 
in  the  turbulent  model.  Because  the  Wilcox  k  -  oo  turbulent  model  (Wilcox,  2006)  is  used  in 
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this  study,  the  conservative  scalars  are  k  and  00.  Because  computational  grids  are  used  to 
represent  the  flow  field,  the  dimension  n  of  the  vector  w  in  the  computation  is  Ngrids  x  Nvars, 
where  Ngrids  is  the  number  of  the  grid  points  and  Nvars  is  the  number  of  flow  variables  at 
each  grid  point.  In  this  study,  Nvars  is  7:  density  (p),  the  three  components  in  momentum 
(pu),  total  energy  [pE),  and  the  two  scalars  in  the  k  -  00  RANS  model  (pk  and  poo).  Similarly, 
the  dimension  of  the  vector  function  F  is  equal  to  n.  Therefore,  the  phase  space  of  this 
dynamical  system  comprises  total  (n  +  1) -dimensional  space.  The  trajectory  of  the  solution 
evolves  in  time  and  in  the  phase  space  from  an  initial  solution  in  w  and  cf>.  In  general,  the 
system  input  parameter  <j)  is  assumed  to  be  given  and  fixed,  and  in  this  case,  an  n- 
dimensional  phase  space  only  for  w  is  considered. 

In  this  study,  we  are  particularly  interested  in  equilibrium  solutions  and  in  the  dynamics 
near  the  equilibria.  An  equilibrium  solution  wo(4>)  at  a  given  parameter  <J>  satisfies 


F(w0;</>)  -  0. 


In  general,  the  equilibrium  solution  set  forms  a  continuous  curve  in  the  phase  space.  The 
dynamics  near  an  equilibrium  solution  (wo;  4>o)  can  be  expressed  in  the  following 
linearized  system: 


dw' 

dt 


(tu0;^o) 


In  this  linearized  system,  w'  £  Rn  is  a  perturbation  vector  from  an  equilibrium  w  0[(J>o), 

where  4>o  is  assumed  to  be  a  constant  value.  *-4(Wo;(|>o)  e  Rnxn  is  the  Jacobian  matrix  of  F, 
evaluated  at  the  equilibrium  solution,  (wo;  c|>o): 


*4(tuo;<!>o)  —  VwF(w0:<f>0) 


*>3  z 


OF  j 
dw , 


(ut0:^o)/ 


l3 


The  eigen-decomposition  of  .4  [wO;4>o]  plays  an  important  role  in  explaining  the  system 

dynamics.  Let  XiM,-  •  »,An  be  eigenvalues  of  [w0;f(>03  and  yi,y2,t  •  .,yn  be  the  associated 
direct  global  modes  [or  right  eigenvectors).  In  other  words,  each  pair  of  Ai  and  yi  satisfies 


,  71. 


Here,  Afs  are  scalars,  and  y/s  are  n-dimensional  vectors.  For  complex  Afs,  the 
corresponding  y/s  are  also  complex.  Without  loss  of  generality,  it  can  be  assumed  that 
Re(Ai)  >  Re[Az)  >•  *  .>  Re[An).  The  real  part  of  an  eigenvalue  stands  for  the  growth  rate  of 
the  corresponding  mode  in  time,  whereas  the  imaginary  part  is  related  to  the  oscillatory 
dynamics.  Therefore,  Ai  is  the  least  stable  because  its  associated  global  mode  grows  the 
fastest  in  time  [if  Re[Ai)  >  0)  or  decays  the  slowest  in  time  [if  Re[Ai)  <  0).  If  all  the 
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eigenvalues  of«'H(Wo;<i>o)  have  a  negative  real  part,  any  point  in  a  neighbor  of  (wo;  4>o)  is 
stable.  On  the  other  hand,  if  few  eigenvalues  have  a  zero  real  part,  the  system  is  marginally 
stable. 

3.3  Methodology  and  configurations 

3.3.1  Methodology  and  configurations 

In  this  study,  the  pseudo-arclength  continuation  method  (Keller,  1977;  Chan  & 
Keller,  1982)  is  applied  to  obtain  the  solution  curve.  In  the  pseudo-arclength  continuation 
algorithm,  an  additional  equation  called  the  “tangential  condition"  is  added  to  F(w;  4>).  To 
obtain  the  (k  +  l)th  solution  point,  two  previous  solutions,  zk  l  =  (wk-1;  <f>kl)  and  zk  =  (wk; 
4>k)  are  first  stored.  Because  the  dimension  of  w  and  <$>  are  n  and  1,  respectively,  the 
dimension  of  z  is  equal  to  n  +  1.  Using  these  two  solution  vectors,  a  tangential  vector  t  =  zk  - 
zk  l  and  the  pseudo-arclength  As  =  ||zk  -  zkl||2  are  calculated.  Here,  ||v||2  stands  for  the  two- 
norm  of  a  vector  v.  Then,  the  next  solution  zk+1  is  searched  on  a  plane  that  satisfies  both  of 
the  following  conditions:  (1)  the  extrapolated  point  zk  +  t  should  be  on  the  plane,  and  (2) 
the  tangential  vector  of  the  plane  is  equal  to  t/As.  In  order  for  zk+1  to  be  on  such  a  plane,  zk+1 
must  satisfy  the  tangential  condition,  T(zk+1,t,  As)  =  tT(z  -zkyAs  -  As  =  0. 

During  the  search  process  for  zk+1,  the  Newton  iterative  method  is  used  (Allgower  & 
Georg,  1997).  When  the  tangential  condition  is  included,  the  non-linear  system  is  expressed 
by  the  vector  function 


F(z) 


G(z)  =l  V-'fA.)  Je  Rn+1, 


and  zk+1  is  the  root  of  G(z)  =  0.  At  each  Newton  iteration  step,  the  solution  is  updated  by  z  := 
z  +  aAz,  where  a  is  a  step  size,  and  Az  is  a  search  direction.  The  range  of  a  is  0  <  a  <  1,  and  it 
is  determined  based  on  the  three-point  parabolic  method  (Kelley,  1995).  The  Az  =  -$z-1G(z), 
where 


R(n+l)x(n+l)^ 


and  ✓Hz  is  the  Jacobian  matrix  of  the  vector  function  F  evaluated  at  z.  The  stopping 
criterion  of  the  Newton  iteration  is  1 1 G (z)  1 1 2  ^  Gabs  or  1 1 G (z)  1 1 2^  |  G (zk  +  t)  1 1 2  <  erei,  where  eabs 
and  €reai  are  the  absolute  threshold  and  the  relative  threshold,  respectively.  In  this  study, 
the  threshold  values  are  eabs  =  10_6  and  erei  =  3  x  1011.  After  satisfying  the  stopping 
criterion,  three  more  Newton  iterations  are  performed  to  ensure  convergence. 

To  prevent  negative  density  and  negative  pressure,  two  barrier  functions,  bi(p)  =  -(Bi  log  p 
and  b2(p)  =  -P2  log  p,  are  added  to  the  mass-conservation  equation  and  the  total-energy- 
conservation  equation  as  source  terms.  The  coefficients  in  the  barrier  functions  are  chosen 
to  be  Pi  =  10'9  and  P2  =  1012  in  this  study,  and  they  become  zero  after  the  stopping 
criterion  of  the  Newton  iteration  is  satisfied.  Because  the  coefficients  are  very  small,  the 
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barrier  functions  are  negligible  in  a  regular  situation.  However,  if  density  or  pressure  nears 
zero  during  the  Newton  iteration,  the  corresponding  barrier  function  gives  a  significantly 
large  function  value,  resulting  in  a  high  slope,  and  thus  the  search  direction  Az  goes  away 
from  the  corresponding  boundary,  p  =  0  or  p  =  0.  The  barrier  functions  with  the  parabolic 
step-size  algorithm  helps  the  Newton  system  become  stable. 

The  flow  solver  that  provides  us  with  F(w;  <$>)  during  the  Newton  iteration  is  our  in-house 
Reynolds-averaged  Navier-Stokes  (RANS)  solver  called  “JOE”  (Pecnik  etal.,  2012).  The  JOE 
solver  is  a  finite-volume  collocated  compressible-flow  solver  using  the  HLLC  shock¬ 
capturing  scheme  (Toro  etal.,  1994)  and  can  handle  unstructured  meshes  in  three- 
dimensional  space.  Its  accuracy  is  second  order  in  space  on  unstructured  meshes  if  shocks 
do  not  occur,  but  the  spatial  accuracy  reduces  to  about  first  order  near  a  shock  due  to  extra 
dissipation  from  the  HLLC  shock-capturing  scheme.  Turbulent  viscosity  in  the  flow  is 
provided  by  the  Wilcox  k  -  co  RANS  model  (Wilcox,  2006). 

The  implementation  of  the  Newton  iteration  shows  a  reasonable  convergence.  The  iterative 

process  converges  in  about  10  iterations  if  the  Jacobian  matrix  *4Z  is  stable,  whereas  more 
than  100  iterations  are  required  if  the  Jacobian  matrix  is  unstable  or  marginally  stable.  This 

is  because  the  condition  number  of  r4  z  grows  as  the  least-stable  eigenvalue  of  the  Jacobian 
matrix  approaches  zero.  When  ,4 z  is  only  marginally  stable,  the  condition  number  of^z 
also  becomes  large,  deteriorating  the  convergence  of  the  Newton  iteration. 

3.3.2  Calculation  of  the  Jacobian  matrix  and  evaluation  of  its  eigen-pairs 

The  Newton  iteration  requires  the  evaluation  of  the  Jacobian  matrix » 4  (w;<w  as  well  as  the 
evaluation  of  the  vector  function  F(w;  4>).  To  compute  *4(W;<j>),  taking  into  account  complex 
geometry  using  unstructured  meshes  as  well  as  shock-capturing  via  the  HLLC  scheme,  we 
employ  the  technique  of  automatic  differentiation  (AD)  (Griewank,  2000).  AD  is  a 
technique  whereby  exact  derivatives  of  a  function  are  calculated  by  computers  without 
truncation  errors,  and  thus  it  is  much  more  accurate  than  traditional  methods  such  as  finite 
differences.  This  technique  has  already  been  used  in  some  fluid  dynamics  applications.  For 
example,  Wang  et  al.  (2012)  applied  AD  to  the  JOE  flow  solver  in  estimating  the  probability 
of  unstart  of  an  inviscid  scramjet  engine.  In  their  study,  AD  played  an  important  role  in 
producing  adjoints  that  were  used  to  reduce  sampling  costs. 

We  rewrote  the  original  code  by  Wang  et  al.  to  obtain  a  better  performance  since  the 
number  of  Jacobian  evaluations  is  several  orders  of  magnitude  higher  than  that  of  Wang 
etal.  (2012)  due  to  the  Newton  iterative  procedures.  To  improve  performance,  the  new 
implementation  takes  advantage  of  modified  data  structures  as  well  as  efficient  access  to 
the  AD  package.  As  a  result,  the  new  code  runs  over  100  times  faster  than  the  original  code 
by  Wang  et  al.  while  giving  the  same  Jacobians. 

The  eigen-decomposition  of  the  Jacobian  matrix  ,4  (w;*)  is  obtained  by  using  the  Arnoldi 
iteration  developed  in  the  SLEPc  parallel  linear-algebra  package  (Hernandez  etal.,  2005), 

because  * 4  [w;(W  is  sparse  and  unsymmetric.  However,  the  convergence  of  the  nominal 
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Arnoldi  method  is  undesirable  when  the  eigenvalues  ofvH(W;<M  are  distributed  in  an  ill- 
favored  way  or  the  size  of  the  Jacobian  matrix  is  not  small  (in  this  study,  the  size  is  about 
1.5  x  106  by  1.5  x  106).  Therefore,  the  shift-and-invert  spectral-transformation 
algorithmQia  &  Zhang,  2002)  is  adopted  for  improved  convergence.  In  this  algorithm,  the 
original  eigen  problem  of  ,4y  =  Ay  is  modified  to 

(-4  -  aiy  'y  =  (A  -  a)~ly. 

The  global-mode  y  in  this  modified  problem  is  the  same  as  the  global  mode  of  the  original 
problem,  and  a  proper  choice  of  the  shift  a  accelerates  the  convergence  of  the  iterative 
process  for  eigenvalues  near  a.  Hence,  in  this  study,  the  least-stable  eigenvalue  Ai  is  first 
found  using  the  nominal  Arnoldi  iteration,  and  other  eigenvalues  are  found  by  the  shift- 

and-invert  spectral-transformation  algorithm  with  a  =  1/2  Ai.  The  matrix  inverse  of  (A  -CTl) 
is  found  by  the  LU  decomposition  that  is  provided  by  the  parallel  direct  sparse  solver, 
MUMPS  (Amestoy  et  al.,  2000). 

3.3.3  Scramjet  configuration 

The  scramjet  model  in  the  study  is  adopted  from  the  model  scramjet  in  the  series  of  the 
HyShot  II  experiments.  The  successful  flight  test  was  performed  at  the  University  of 
Queensland  (Smart  etal.,  2006)  followed  by  ground  experiments  both  at  the  University  of 
Queensland  (Frost  etal.,  2009)  and  in  the  High  Enthalpy  shock  tunnel  Gottingen  (HEG) 
facility  of  the  German  Aerospace  Center  (DLR)  (Gardner  et  al.,  2004;  Hannemann 
et  al.,  2010;  Laurence  et  al.,  2013,  2015).  In  the  HyShot  II  flight  experiment,  the  free-stream 
properties  are  M  =  7.8,  T  =  242  K,  and  p  =  1711  Pa  when  HyShot  II  was  in  an  altitude  of  27 
km,  and  the  angle  of  attack  of  the  vehicle  at  this  location  was  3.6°.  These  free-stream 
conditions  were  also  investigated  in  Gardner  etal.  (2004),  and  the  later  experimental 
studies  at  DLR  aimed  at  an  altitude  of  28  km,  using  similar  free-stream  conditions.  This 
series  of  experimental  studies  revealed  data  sets  in  unstart  physics,  but  the  detailed 
analysis  of  the  scramjet  model  still  relied  on  ID  analysis  or  simulation  data  of  the  same 
geometry  because  access  to  the  comprehensive  flow  field  was  limited  in  the  experiments. 

The  scramjet  model  in  the  experiments  includes  an  18°  intake  ramp,  a  0.3  m-long  constant- 
area  combustor,  and  an  expansion  nozzle.  However,  because  the  supersonic  flow  on  the 
intake  ramp  is  not  affected  by  the  downstream  changes  in  the  combustor-nozzle  system 
during  the  unstrart  process,  the  intake  ramp  is  not  included  in  this  study.  Figure  20  gives 
the  comparison  between  the  full  geometry  and  the  reduced  geometry.  The  contour  levels  in 
the  figure  represent  velocity  divergence,  showing  compression  regions  as  black  and 
expansion  regions  as  white.  Therefore,  shock  waves  are  shown  as  black  lines,  whereas 
white  regions  show  expansion  waves.  Figure  20(a)  shows  the  two-dimensional  flow  field  of 
the  full  geometry  with  (J>  =  0,  and  the  computational  mesh  is  taken  from  the  computational 
study  by  Pecnik  et  al.  (2012)  that  considered  the  entire  intake-combustor-nozzle  system.  In 
this  full-geometry  calculation,  the  properties  of  the  free  stream  at  the  domain  inlet  are 
matched  with  the  corresponding  values  in  the  ground  experiment  by  Gardner  et  al.  (2004), 
and  the  walls  are  isothermal  at  300  K.  The  intake  ramp  starts  at  x  =  0  m,  where  x  is  the 
horizontal  coordinate,  and  the  leading  edge  of  the  lower  wall  of  the  combustor  is  located  at 
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Xiniet  =  0.35  m.  Because  the  combustor  is  0.3  m  long,  the  following  expansion  nozzle  starts  at 
x throat  =  0.65  m.  On  the  other  hand,  the  inlet  of  this  reduced  domain  is  located  at  about  4.3 
mm  upstream  of  Xiniet,  as  shown  in  Figure  20(b).  On  the  combustor  walls,  the  grid  size  in  the 
wall-normal  direction  is  fine  enough  to  capture  the  growth  of  boundary  layers  on  these 
walls.  The  number  of  grid  points  in  this  reduced  domain  is  0.21  x  106. 


Fig.  20.  Divergence  contours  at  4>  =  0:  (a)  full  computational  domain;  (b)  reduced  domain. 
Black  regions  indicate  negative  divergence  (compression);  white  regions  indicate  positive 
divergence  (expansion).  The  scales  in  (a)  and  (b)  are  not  the  same. 


The  flow  profiles  at  the  domain  inlet  in  Figure  21(b)  are  taken  from  the  same  grid  locations 
in  the  full  domain  simulation  in  Figure  21(a),  and  the  flow  profiles  are  applied  to  the 
reduced  domain  as  the  Dirichlet  boundary  condition.  Figure  21  highlights  the  two 
computational  domains  near  Xiniet-  The  Mach  number  at  the  reduced-domain  inlet  is 
approximately  M  =  2.62,  and  the  flow  direction  becomes  parallel  to  the  upper  wall  and 
lower  wall  of  the  combustor  before  entering  into  the  combustor.  At  the  inlet  of  the 
combustor  (xiniet  =  0.35  m),  two  oblique  shocks  are  formed  due  to  the  round  shape  of  the 
lower-wall  leading  edge.  The  upper  oblique  shock  that  comes  into  the  combustor  reflects 
from  the  combustor  walls,  generating  an  oblique-shock  train. 


(a)  (b) 


0.35  0.4 


Fig.  21.  Divergence  contours  near  the  combustor  inlet  at  (j)  =  0:  (a)  full  computational 
domain,  (b)  reduced  domain.  Black  regions  indicate  negative  divergence  (compression); 
white  regions  indicate  positive  divergence  (expansion). 


In  this  study,  the  heat  from  combustion  is  released  to  the  scramjet  system  by  a  surrogate 
model.  In  the  HyShot  II  experiments,  hydrogen  fuel  was  injected  in  the  wall-normal 
direction  at  58  mm  downstream  of  the  leading  edge  of  the  combustor,  and  the  fuel  was 
mixed  with  the  incoming  air  and  combusted  afterward.  The  heat  release  from  combustion 
is  followed  by  a  pressure  rise  along  the  stream-wise  direction,  and  excessive  heat  release 
results  in  unstart.  However,  capturing  detailed  combustion  chemistry  requires  excessive 
computational  costs,  and  thus  a  low-order  model  is  required  when  a  large  number  of 
computations  is  needed  to  obtain  the  solution  curve.  Low-order  models  have  been 
proposed  in  scramjet  studies,  and  many  of  them  are  one-dimensional  models.  For  example, 
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Mitani  etal.  [2003]  used  one-dimensional  simulations  to  predict  thrust  of  two  different 
scramjets,  and  the  thrust  prediction  by  the  simple  simulations  agreed  reasonably  well  with 
the  experimental  measurements  until  the  fuel  equivalence  ratio  exceeded  the  unstart  limit. 
In  the  study  of  HyShot  II  by  Laurence  et  al.  (2013),  the  unstart  limit  in  the  equivalence  ratio 
was  successfully  predicted  by  a  one-dimensional  model  based  on  Rayleigh’s  flow 
assumptions,  and  the  authors  reported  that  the  onset  of  unstart  is  sensitive  to  the  total  heat 
release  from  the  combustion.  Doolan  &  Boyce  (2008)  used  a  quasi  one-dimensional  mixing- 
and-combustion  model  to  estimate  the  performance  of  the  ground  experiment  by  Boyce 
etal.  (2000),  and  the  model  output  agreed  well  with  the  experimental  data.  Similarly, 
Tourani  (2011)  showed  that  one-dimensional  simulations  can  closely  predict  the  overall 
evolution  of  the  flow  in  the  scramjet  combustor  as  shown  in  the  study  of 
Oevermann  (2000). 

The  heat-release  model  in  this  study  is  the  model  proposed  by  Wang  etal.  (2012).  In  this 
model,  heat  is  added  to  the  system  through  a  volumetric  term  given  in  the  following 
formula: 


Q  =  <PfstHfm<lirr)(x/Lc ), 

where  4>  is  the  heat-release  rate,  and  q(xlc)  =  1  -  e  tCcxlclDc  is  the  heat  distribution  function 
along  the  stream-wise  direction.  The  parameters  in  the  above  equation  are  given  in  Table  2. 
Using  this  model,  the  heat  release  from  the  combustion  in  the  HyShot  II  experiment  is 
modeled,  whereas  the  boundary  layer  and  the  oblique  shock  train  in  the  combustor  are 
captured  in  detail  in  the  two-dimensional  flow  simulations. 


mbol 

Definition 

Baseline  value 

f. 

Stoichiometric  fuel/air  ratio 

0.028 

Hr 

Fuel  heating  value  of  hydrogen  fuel 

120  MJ/kg 

L 

Combustor  length 

0.368  m 

K 

Fraction  of  completed  combustion 

0.95 

Dc 

Shape  parameter 

0.75 

a 

Shape  parameter 

-log{  1  -  X),» 

X 

Combustion  ignition  position 

0.418  m 

Table  2.  Parameters  used  in  the  heat-release  model 


3.4  Results 

3.4.1  Solution  curve 

The  solution  curve  of  the  two-dimensional  HyShot  II  is  calculated  using  the  pseudo- 
arclength  continuation  technique,  and  a  part  of  the  curve  is  given  in  Figure  22.  The 
horizontal  axis  of  the  plot  is  the  heat-release  rate  4>  given  in  Eq.  (7),  and  the  vertical  axis  is 
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the  average  density  in  the  combustor-nozzle  system  that  was  an  effective  metric  to 
represent  the  unstart  in  a  simple  converging-diverging  nozzle  in  Jang  etal.  [2012).  The 
numerical  continuation  starts  from  a  converged  steady-state  RANS  solution  at  cj>  =  0.3970 
using  the  nominal  JOE  solver,  and  the  solution  curve  proceeds  with  the  aid  of  the  numerical 
continuation  method.  The  initial  step  size  in  <J>  is  A<J>  =  2.382  x  10  3,  and  Acf)  changes 
adaptively  at  each  point  on  the  solution  curve  to  get  the  minimum  number  of  iterations  in 
the  Newton  method.  Near  P2  in  Figure  22,  for  example,  A(J>  is  reduced  to  8.922  x  10  5 
because  the  Newton  method  requires  an  excessive  number  of  iterations  if  A<J>  is  not 
reduced  near  this  solution  point. 
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Fig.  22.  Solution  curve  of  the  300  mm  case. 

When  (J>  is  low,  the  solution  curve  in  Figure  22  is  almost  a  linear  function  of  (J>.  However,  as 
the  solution  proceeds  on  the  solution  curve  from  PI  to  P2  in  Figure  22,  the  slope  of  the 
curve  becomes  steeper,  forming  an  inflection  point  at  P2. 


At  the  infection  point  P2,  a  weak  shock  structure  is  first  formed  near  the  throat,  and  as  cf) 
increases  from  P2,  the  shock  structure  moves  upstream  and  becomes  stronger.  The 
divergence  contours  near  the  throat  (xth)  at  the  four  solution  points  shown  in  Figure  23 
show  the  evolution  of  the  shock  structure.  The  enlarged  images  in  Figure  23  clearly  show 
the  evolution  of  the  normal  shock  near  the  throat.  In  the  divergence  contours  at  P2  given  in 
Figure  23(b),  a  weak  shock  structure  is  formed  near  x  =  0.643  m,  that  is  not  shown  at  a 
lower  4>  (e.g.,  the  divergence  contours  at  PI  given  in  Figure  23(a)).  If  a  higher  value  of  (J>  is 
applied,  the  normal  shock  becomes  stronger,  and  its  position  moves  upstream,  as  shown  in 
Figure  23(c)  and  (d). 
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This  result  is  similar  to  the  observations  in  the  ground  experiments  by  Laurence 
etal.  [2013).  In  their  experiment,  pressure  was  measured  by  pressure  tabs  on  both  the 
injector-side  wall  and  the  cowl-side  wall,  and  the  location  of  a  sudden  rise  in  the  wall- 
pressure  distribution  was  also  stationary  in  time.  The  sudden  rise  in  wall-pressure  is 
believed  to  be  related  to  the  leading  shock  in  the  unstart  structure,  and  the  location  of  the 
sudden  rise  was  a  function  of  the  fuel  equivalence  ratio  in  the  experiment  (4>exP).  When  c()exp 
^  0.66,  for  example,  the  location  of  the  leading  shock  moved  to  just  downstream  of  the 
injector,  and  there  was  a  noticeable  uncertainty  of  the  leading-shock  location  at  a  high 
value  of  (j>exP  in  their  experiments.  If  c|>exp  increases  to  about  1.1,  the  location  of  the  pressure 
jump  passes  the  injector,  but  the  exact  position  was  not  reported  because  pressure  in  the 
upstream  from  the  injector  was  not  measured. 


[a)  PI  fd>  =  0.5497)  (b)  P2  [<b  =  0.5611) 


063  064  065  066  063  064  065  066 


(c)  P3  [4>  =  0.5699)  [d)  P4  [cj>  =  0.5733) 

■ 

0.63  0.64  0.65  0.66  0.63  0.64  0.65  0.66 

Fig.  23.  Velocity  divergence  contours  near  xth  =  0.65:  (a)  PI  in  Figure  22  (4>  =  0.5497),  (b) 
P2  [<£  =  0.5611),  [c)  P3  [4>  =  0.5699),  (d)  P4  (4>  =  0.5733).  Black  regions  indicate  negative 
divergence  [compression);  white  regions  indicate  positive  divergence  [expansion). 

Figure  24  shows  the  Mach  number  contours  at  the  four  different  points  on  the  solution 
curve.  In  the  Mach  number  contours  at  PI,  the  flow  is  supersonic  in  most  of  the  region 
between  the  walls,  and  the  subsonic  zones  are  confined  in  the  boundary  layer  on  the  walls. 
When  4>  is  increased  to  0.5611  [P2),  the  subsonic  zone  becomes  thicker,  especially 
downstream  of  the  normal  shock.  If  more  heat  is  added  to  the  system,  the  subsonic  zones 
grow  to  the  centerline  of  the  combustor,  and  they  finally  form  a  full  subsonic  band  that 
spans  the  total  height  of  the  duct.  The  first  formation  of  the  full  subsonic  band  is  found  at 
P3,  and  a  wider  band  is  observed  at  a  more  upstream  location  as  4>  increases. 
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Fig.  24.  Mach  number  contours  near  xth  =  0.65  at  four  different  points  on  the  solution  curve 
given  in  Fig  3:  [a)  PI  in  Figure  22  =  0.5497),  (b)  P2  (cf>  =  0.5611),  (c)  P3  ((j)  =  0.5699),  (d) 

P4  (cf)  =  0.5733).  White  regions  are  subsonic  (M  <  1). 

3.4.2  Linear  dynamics 

The  eigenvalues  of  the  Jacobian  matrix, 4 at  four  points  (PI,  P2,  P3,  and  P4)  on  the 
solution  curve  are  given  in  Figure  7.  Because  eigen-decomposition  requires  significant 
computing  resources,  only  13  eigenvalues  with  the  largest  real  part  (Ai,,  •  ,,Ai3)  are 
calculated.  Of  the  13  eigenvalues,  5  eigenvalues  are  purely  real,  whereas  the  other  8 
eigenvalues  form  four  complex-conjugate  pairs.  At  PI,  the  real  part  of  the  least-stable 
eigenvalue  Ai  is  about  -1.0  x  10'6,  whose  absolute  value  is  an  order  of  magnitude  higher 
than  the  Ai’s  at  the  other  three  solution  points.  In  addition,  Ai  at  PI  is  close  to  the  other 
eigenvalues  (A2,*  •  -,Ai3),  but  when  <J>  increases  to  the  value  at  P2,  Ai  is  separated  from  the 
other  eigenvalues,  approaching  zero.  At  P3  and  P4,  the  distance  between  Ai  and  A2  is  still 
considerably  larger  than  that  at  PI.  Therefore,  as  cf)  increases,  the  mode  associated  with  Ai 
can  survive  longer,  while  the  other  modes  decay  quickly.  This  means  that  this  slow 
dynamics  of  Ai  is  associated  with  the  evolution  of  the  moving-shock  structure. 
Furthermore,  Ai  at  P2  comes  close  to  zero,  and  thus  the  system  is  marginally  stable  when 
the  shock  structure  is  first  formed. 


The  least-stable  direct  global  mode  at  the  inflection  point  (P2)  is  shown  in  Figure  26  along 
with  the  divergence  contours  in  Figure  26(a)  that  indicate  the  location  of  the  normal  shock 
at  x  ^  0.642.  In  the  contour  plots  in  the  Figure  26(b),  (c),  and  (d),  the  amplitude  and  the 
sign  are  arbitrary  because  of  normalization.  Figure  26(b)  exhibits  the  vector  field  of 
(yi,ui.yi, 112)  =  (yi,Pui/p,yi,pu/p).  The  vector  field  shows  two  big  circular  motions  behind  the 
shock,  and  they  are  symmetric  with  respect  to  the  centerline  between  the  upper  and  lower 
walls.  The  circular  motions  near  the  walls  are  directed  to  downstream,  and  they  merge  into 
one  upstream  motion  on  the  centerline.  Figure  26(c)  and  (d)  show  the  contour  plots  of  yi,p 
and  yi.pE,  respectively.  In  both  plots,  a  noticeable  peak  is  found  on  the  centerline 
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downstream  of  the  normal  shock,  and  the  shape  of  the  peak  is  an  oval  whose  long  axis  is 
oriented  along  the  wall-normal  direction. 


10x10_7: 

D  * 

□  X 

5 

\ 

!  ° 

-5 

□  X 

-10x1  O'7 

□  * 

-2x1  O'6  -1  0 


Re(A) 


O  PI 
A  P2 
□  P3 
x  P4 


Fig.  25.  13  eigenvalues  with  largest  real  parts  at  4  different  points  on  the  solution  curve  of 
Fig  3.  Black  circles  =  PI;  red  triangles  =  P2;  green  squares  =  P3;  orange  crosses  =  P4. 
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Fig.  26.  The  least-stable  global  mode  at  the  inflection  point  P2  (enlarged  near  the  throat): 
(a)  divergence  of  velocity  field  in  the  base  flow  field  (black=compression, 
white=expansion),  (b)  vector  field  of  (yi,ui,yi,u2),  (c)  density  yi,p,  (d)  total  energy  yi,PE. 
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3.4  Conclusions 


The  equilibrium  solution  curve  of  the  two-dimensional  HyShot  II  is  found  as  a  function  of 
the  heat-release  rate  4>  by  the  pseudo-arclength  technique  with  Newton  iteration.  At  the 
inflection  point  on  the  solution  curve,  a  shock  structure  is  first  found,  and  the  shock 
structure  moves  upstream  with  increasing  <$>,  as  previously  observed  in  a  ground 
experiment  of  HyShot  II.  The  linear-system  analysis  reveals  the  separation  of  slow  and  fast 
dynamics,  and  visualization  of  the  least-stable  global  mode  at  the  inflection  point  depicts 
the  strengthening  or  weakening  mechanisms  of  the  shock  structure. 

Based  on  the  onset  dynamics  of  the  shock  structure  found  in  this  study,  unstart-mitigation 
mechanisms  can  be  suggested  in  connection  with  the  linearized  dynamics.  In  particular,  the 
adjoint  global  modes  of  the  linearized  system  are  related  to  the  receptivity  of  the  system  to 
perturbations.  Therefore,  the  optimal  control  mechanisms  can  be  found  by  considering 
both  the  direct  global  modes  and  the  adjoint  global  modes.  Calculating  the  adjoint  global 
modes  and  finding  the  unstart-mitigation  mechanisms  are  the  subjects  of  the  ongoing 
study. 
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